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Abstract.  Gravity  waves  on  the  sea  surface  reflect  coupling  between  the  atmo¬ 
sphere  and  ocean.  They  play  important  dynamical  roles  in  the  transfer  of  momen¬ 
tum.  heat,  and  scalars,  but  how  they  interact  with  background  oceanic  (and  atmo¬ 
spheric)  turbulence  is  still  not  well  understood.  This  is  especially  true  under  high 
winds  where  waves  foster  Langmuir  circulations  but  also  break  intermittently. 
Results  from  high-resolution  large  eddy  simulations  (LESs)  that  include  stochas¬ 
tic  breaking  waves  and  the  Craik-Leibovich  “vortex  force”  are  analyzed.  Under 
certain  sea  states  turbulence  (anisotropic  seed  vorticity)  generated  by  breaking 
waves  initiates  Craik-Leibovich  instabilities  that  lead  to  coherent  structures,  viz., 
Langmuir  cells  and  intermittent  downwelling  jets  both  of  which  impact  ocean 
mixing.  We  suggest  a  mechanism  as  to  how  breaking  waves  couple  with  the 
vortex  force  that  leads  to  the  formation  of  these  flow  structures. 

1.  Introduction 

Gravity  waves  at  the  sea  surface  impact  momentum  and  scalar  exchange  in 
the  marine  boundary  layers.  Waves  are  particularly  important  in  the  upper  ocean 
boundary  layer  (OBL)  as  they  break  intermittently  and  thus  sporadically  inject 
momentum  into  the  water  column  leading  to  the  growth  of  currents  [1;  2],  This 
momentum  (and  energy)  transfer  from  winds  to  waves  to  currents  is  beautifully 
captured  in  photographs  taken  from  aircraft  under  high  wind  conditions  (figure  1). 
Also,  when  waves  and  currents  are  averaged  over  multiple  periods,  a  net  wave- 
current  interaction  fosters  the  development  of  Langmuir  circulations  [3;  4],  The 
elegant  and  delicate  mathematical  steps  and  theoretical  arguments  leading  to  these 
conservative  wave-current  interactions  appear  as  “vortex  forces”  in  the  equations 
of  motion  [5;  3;  6;  7].  There  has  long  been  intuitive  speculation  that  wave  breaking 
and  Langmuir  circulations  interact  [8;  9;  4],  but  observing  and  elucidating  this 
coupling  is  difficult  in  complex  sea  states  owing  to  the  intermittent  nature  and 
widely  varying  scales  associated  with  breaking  waves.  Also,  Stokes  drift  which 
contributes  to  the  vortex  force  is  a  Lagrangian  property  of  the  wavefield  and  not 
readily  deduced  from  fixed-point  Eulerian  measurements. 

Recently  we  have  developed  idealized  large-eddy  simulation  (LES)  models  of 
the  turbulent  OBL  which  include  both  phase-averaged  wave-current  interactions 
and  intermittent  wave  breaking  [10;  11;  12;  13].  These  turbulence  resolving  mod¬ 
els  can  be  used  to  examine  possible  couplings  between  wave  breaking  and  Lang¬ 
muir  circulations.  In  our  LES  model,  the  wave-averaged  dynamical  effects  on 
currents  are  the  result  of  an  asymptotic  theory  pioneered  by  Craik  and  Leibovich 
[5;  3]  and  later  extended  by  McWilliams  et  al.  [6;  7].  Effects  of  intermittent 
breaking  are  included  as  a  probability  density  function  (PDF)  of  parameterized 
stochastic  impulses  that  depend  on  time  and  space.  The  PDF  of  breaking  varies 
with  wind  speed  and  wave  age  and  obeys  the  conservation  rule  that  when  aver¬ 
aged  over  long  time  and  large  area  reduces  to  a  traditional  bulk  formula  for  con¬ 
stant  wind  stress  T.  Wave  effects  appear  in  our  LES  model  equations  for  currents. 
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Figure  1 :  Photograph  of  breaking  waves  in  the  Gulf  of  Tehuantepec  at  a  wind 
speed  of  approximately  20ms  1 .  The  image  is  from  a  low-level  flight  of  the 
NCAR  C-130.  Courtesy  W.  Melville  and  C.  Friehe. 


scalars,  and  subgrid-scale  turbulent  kinetic  energy  (TKE).  Schematically,  the  LES 
equation  for  resolved-scale  momentum  p u  with  wave  effects  is 

37  =  ...  +  u5'  xro  +  £  A(,)  ,  (1) 

dr 

where  spatially  filtered  variables  are  indicated  by  ( ).  In  (1),  the  resolved  velocity 
is  m,  the  vorticity  is  ©  =  V  x  m,  us  is  the  Stokes  drift  associated  with  the  wave 
field,  and  A*')  are  compact  (in  time  and  space)  accelerations  designed  to  mimic 
wave  breaking.  At  any  time  step  of  the  simulation  the  number  of  active  impulses 
n  varies  depending  on  the  PDF  of  breaking  waves.  Dots  in  (1)  denote  all  other 
terms  in  conventional  LES  of  geophysical  boundary  layers,  viz.,  advection,  buoy¬ 
ancy,  Coriolis,  pressure  gradients,  and  subgrid-scale  flux  divergence,  for  example 
see  [14],  The  vertical  profile  of  Stokes  drift  is  a  primary  input  to  LES  of  the  OBL 
with  wave-averaged  effects  when  currents  are  generated  by  uniform  stress.  Typ¬ 
ically  uSt(z)  is  estimated  from  a  single  monochromatic  wavenumber  component 
k  selected  to  match  the  peak  in  the  wave  height  spectrum;  then  the  Stokes  drift 
has  depth  variation  uSt  ~  e2kz  [13;  15;  16;  17],  In  these  previous  LES  the  vortex 
force  has  been  shown  to  have  a  significant  influence  on  the  turbulent  eddies  and 
their  vertical  fluxes  in  the  OBL  through  the  generation  of  Langmuir  circulations 
[5;  13]. 

Our  current  LES  model  of  the  OBL  builds  on  past  work  and  incorporates  ad¬ 
ditional  wave  effects.  The  most  important  are:  (1)  the  vertical  profile  of  Stokes 
drift,  which  appears  in  the  vortex  force,  is  obtained  by  integration  of  empirical 
wave  height  spectra  [18]  over  all  wavenumbers  k\  and,  (2)  the  assumption  of  con¬ 
stant  stress  is  replaced  by  stochastic  forcing.  As  a  result  of  these  and  other  im¬ 
provements,  the  LES  model  accounts  for  varying  sea  state  expressed  in  terms  of  a 
bulk  wave  age  cp/uta  where  cp  is  the  phase  speed  of  the  peak  in  the  wave  height 
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Figure  2:  Flow  visualization  of  the  m—  current  near  the  water  surface  (z  =  —  1 ,0m) 
from  an  LES  with  stochastic  wave  breaking  and  vortex  force  at  a  wind  speed  of 
SOms^1  [10].  The  horizontal  (Ax,  Ay)  resolution  of  this  calculation  (500  x  500 
gridpoints)  is  1.5m.  The  vertical  spacing  (128  gridpoints)  varies  from  1.0m  at  the 
surface  to  1.2m  at  the  bottom  of  the  computational  domain.  The  positive  (red) 
contours  are  locations  of  breaking  waves.  The  color  bar  is  in  units  of  ms-1  and 
the  winds  are  from  left  to  right. 


spectrum  and  uta  is  the  atmospheric  friction  velocity.  A  typical  LES  flow  field 
from  recent  calculations  at  a  wind  speed  Ua  =  30ms  1  approaching  wind-wave 
equilibrium  =  30  is  shown  in  figure  2.  This  snapshot  of  the  u— current 

is  rich  in  detail;  a  spectrum  of  breaking  waves  in  various  stages  of  development 
is  visible,  and  the  range  of  breaker  scales  varies  from  small  breakers  0(1  )m  to 
large  breakers  0(50  —  75)m.  We  note  the  flow  structures  shown  in  figure  2  are 
not  present  in  comparable  simulations  driven  by  constant  stress.  Our  focus  is  on 
analyzing  the  vortex  force  V f  =  u5'  x  to  from  LES  solutions  for  varying  wave  age 
when  x  is  replaced  by  A'1'1,  i. e. ,  in  the  presence  of  intermittent  wave  breaking  as 
depicted  in  figure  2. 


2.  Results 


For  a  wave  field  propagating  solely  in  the  alongwind  direction  (x— direction) 
us  =  uSt(z) i  and  then  the  vortex  force  has  only  two  non-zero  components 

Vf  =  K*(0i,  COyk)  (2) 


where  the  spanwise  and  vertical  vorticity 


C0V 


and  to,. 


/  dv  du\ 
\dx  dy  J 


(3) 


With  an  externally  imposed  Stokes  profile  the  vortex  force  is  simply  the  resolved 
vorticity  field  weighted  by  an  exponentially  decaying  function  of  z.  In  our  model 
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Figure  3:  Vertical  profiles  of  root-mean-square  (rms)  values  of  the  vortex  force  in 
the  presence  of  breaking  waves  from  LES  at  a  wind  speed  Ua  =  15  ms  *.  Panels 
a)  and  b)  are  (' Vfy,Vfz )  components,  respectively.  Wave  age  cp/u*a  =  [19,23,30] 
(red,  blue,  green)  curves  respectively.  No  breaking  is  solid  black  curve.  The  vortex 
force  is  normalized  by  the  water  friction  velocity  and  mixed  layer  depth  ula/\h\. 


there  is  no  feedback  mechanism  allowing  breakers  to  modify  Stokes  drift,  but 
breaking  does  interact  with  the  vortex  force;  i.e.,  in  (1)  each  breaker  impulse  A'" 
generates  ©  capable  of  coupling  to  the  vortex  force  as  indicated  by  (2).  For  winds 
and  waves  approaching  equilibrium  our  PDF  of  breaking  emphasizes  small  scales, 
whereas  for  strongly  forced  growing  seas  is  weighted  more  heavily  towards  the 
peak  in  the  wave  height  spectrum  [10].  The  scale  and  magnitude  of  the  vorticity 
generated  by  breaking  thus  depend  on  wave  age  [19]  and  wind  speed  [20]  which 
potentially  makes  the  vortex  force  intermittent  and  sea  state  dependent. 


Figure  4:  Vertical  profiles  of  the  (y,z)—  components  of  the  vortex  force  flatness 
for  the  same  simulations  as  in  figure  3.  a)  is  the  y— component  and  b)  is  the 
Z— component  of  the  vortex  force,  respectively. 


In  figures  3  and  4  we  show  vertical  profiles  of  the  root-mean-square  (rms)  and 
flatness  (or  kurtosis)  of  V  f.  (Statistics  at  each  z  are  obtained  by  averaging  across 
x  —  y  planes  and  over  20  3D  volumes.)  First  we  notice  that  overall  the  strength 
of  breaking  impacts  the  statistical  moments.  In  all  cases  the  fluctuations  in  vortex 
force  are  large  compared  to  the  average  surface  stress,  hence  both  components 
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Figure  5:  Probability  density  function  of  the  vortex  force  in  the  presence  of  break¬ 
ing  waves  at  z  =  —  1 .33m.  Wave  age  cp/ uta  =  [19, 23, 30]  (red,  blue,  green)  curves 
respectively.  No  breaking  is  solid  black  curve  and  the  dotted  line  is  a  Gaussian 
distribution.  Panels  a)  and  b)  are  for  the  components  (Vfy,Vfz),  respectively. 


of  Vf  are  dynamically  important  compared  to  shear  forces  as  expected  for  Lang¬ 
muir  turbulence  [17;  13].  Formally  the  vertical  component  Vrz  is  analogous  to  a 
buoyancy  force  in  the  equations  of  motion  with  its  mean  value  absorbed  in  the 
hydrostatic  balance  [3],  but  the  fluctuations  in  the  vertical  component  of  vortex 
force  are  important.  Surprisingly,  near  the  water  surface  the  rms  intensity  of  Vjy 
decreases  in  the  presence  of  breaking  while  the  vertical  component  V fz  is  largely 
unaffected.  A  spectral  analysis  of  the  vortex  force  (not  presented)  shows  that  these 
trends  hold  across  spatial  scale.  Previously  we  reported  [10]  that  the  main  impact 
of  breaking  is  to  weaken  the  coherence  of  Langmuir  cells  accompanied  by  small 
reductions  in  the  spanwise  variance  v'2  compared  to  simulations  driven  by  con¬ 
stant  stress  alone.  This  translates  into  lower  values  of  V',y  leaving  VL  unaffected. 
This  is  expected  based  on  the  definitions  (2)  and  (3).  Digging  deeper  into  the 
statistics  we  notice  that  breaking  alters  the  flatness  of  Vf  as  shown  in  figure  4. 
The  large  departures  from  a  Gaussian  distribution  (flatness  values  =  3)  indicate 
that  the  vortex  force  has  taken  on  a  spotty  spatially  intermittent  character  likely 
mimicking  the  intermittency  of  the  breaker  distribution  as  shown  in  figure  2.  The 
intermittent  character  of  the  vortex  force  is  further  quantified  in  the  PDFs  shown  in 
figure  5.  The  PDF  of  Vry  is  symmetrical,  non-Gaussian,  and  exhibits  long  tails  in 
the  presence  of  strong  breaking.  Meanwhile  the  PDF  of  the  vertical  vortex  force 
is  clearly  asymmetrical,  being  biased  towards  positive  values  in  the  presence  of 
vigorous  breaking. 

3.  Discussion 

These  statistical  distributions  highlight  that  breaking  and  vortex  force  can  cou¬ 
ple  in  our  LES.  Our  interpretation  of  how  this  coupling  occurs  is  based  on  labora¬ 
tory  observations  of  breaking  waves  [2]  and  the  Craik-Leibovich  instability  mech¬ 
anism  CL2  [3],  In  figure  6  we  show  the  time  and  space  evolution  of  the  flowfields 
induced  by  an  idealized  2D  breaker  in  a  laboratory  wave  tank.  The  aftermath  of 
this  breaking  event  is  clear.  A  coherent  spanwise  vortex  is  created  by  the  overturn¬ 
ing  breaker,  and  as  time  advances  it  propagates  downstream  and  descends  deeper 
into  the  water.  The  lifetime  of  this  vortex  is  long,  more  than  50  wave  periods. 
Based  on  these  measurements  we  conclude  that  the  vorticity  field  generated  by 
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Figure  6:  Streamlines  illustrating  the  mean  flow  development  in  the  aftermath  of 
a  2D  breaking  wave  from  laboratory  measurements  of  Melville  et  al.  [2],  The 
color  code  shows  the  magnitude  of  the  mean  velocity.  The  main  feature  is  a  co¬ 
herent  vortex  with  axis  aligned  with  the  y—  span wise  direction.  It  slowly  propa¬ 
gates  downstream  and  deepens  as  time  advances.  From  top-to-bottom  the  non- 
dimensional  time  t  =  [3,10.5,26.5,34.5,42.5,50,58].  t  and  streamwise  distance 
x  are  made  dimensionless  by  the  pre-breaking  wave  period  and  wavelength. 


6 


Vortex  generated  by  a  breaking  wave 


Figure  7:  Sketch  illustrating  the  CL2  instability  mechanism  for  breaker-induced 
generation  of  Langmuir  circulations,  adapted  from  [3],  The  initial  forward  im¬ 
pulse  from  a  3D  breaking  wave  generates  a  surface  current  perturbation  u(y,z), 
spanwise  vorticity  Cby  beneath  the  water  surface,  and  vertical  vorticity  ftj-  at  the 
spanwise  edges  of  the  breaker,  ftj-  is  amplified  by  the  vortex  force  Vfy  and  en¬ 
hances  the  CL2  instability.  Positive  signed  (by  generated  by  breaking  biases  the 
vertical  component  Vrz  towards  positive  values. 


breaking  contains  a  coherent  anisotropic  flow  structure.  In  figure  7  the  essential 
ingredients  of  the  CL2  mechanism  are  presented.  We  have  modified  the  origi¬ 
nal  Craik-Leibovich  interpretation  to  include  the  idealized  single  breaking  wave 
depicted  in  figure  6. 

A  scenario  of  the  CL2  instability  process  instigated  by  breaking  is  as  follows: 
Initially  breaking  supplies  a  forward  impulse  u(y,z)  that  leaves  behind  a  vorticity 
field  with  preferred  orientation  and  sign  (covj,(0:k).  Along  the  centerline  of  the 
perturbation  v  =  0,  (0  — >  (by  j  and  the  sign  of  the  vortex  rotation  is  positive  as  in 
figure  6.  However  the  breaking  process  is  taken  to  be  3D  in  our  simulations,  i.e., 
breakers  span  a  finite  length  A  in  the  y— direction.  As  vortex  lines  are  continuous, 
at  the  edges  y  =  ±A/2  the  spanwise  breaker  vortex  must  bend  vertically  upward 
towards  the  water  surface  as  depicted  in  figure  7;  thus  a  3D  breaking  wave  gen¬ 
erates  significant  amounts  of  vertical  vorticity  ±(bz  as  well  as  spanwise  vorticity 
(by.  Given  this  seed  perturbation  the  CL2  mechanism  then  becomes  activated. 

The  breaker  induced  vorticity  components  directly  contribute  to  the  vortex 
force  as  indicated  by  (2).  We  immediately  notice  with  (by  >  0  a  positive  signed 
vertical  vortex  force  V/z  is  generated  along  the  centerline  v  =  0.  The  PDF  distri¬ 
bution  in  figure  5b  reflects  this.  It  is  strongly  asymmetrical  with  a  bias  towards 
positive  values  in  the  presence  of  large  scale  breaking.  Vertically  oriented  vorticity 
is  crucial  to  the  CL2  mechanics.  Since  the  sign  of  the  vertical  vorticity  rotation  is 
opposite  at  the  spanwise  edges  of  each  breaker,  equal  strength  but  opposing  vortex 
forces  Vfy  =  uSt  cb,  are  generated  aty  =  ±A/2.  Our  idealization  of  the  CL2  process 
is  symmetrical  about  the  plane  y  =  0  which  is  reflected  in  the  symmetrical  PDF  of 
Vfy  shown  in  figure  5.  Opposing  vortex  forces  V/y  promote  the  CL2  instability.  As 
time  advances  from  the  initial  breaker  impulse  the  two  opposing  spanwise  vortex 
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Figure  8:  Temporal  evolution  of  breaker  flowfields  under  the  action  of  the  vortex 
force  at  z  =  —1.14m  for  a  simulation  with  wave  age  cp/u*a  =  19.  The  resolved 
vertical  vorticity  CO-  is  shown  in  the  color  images  and  the  resolved  streamwise 
current  u  as  solid  contour  lines.  The  color  bar  at  the  top  of  the  figure  is  in  units 
of  s_1  and  the  contour  lines  are  evenly  spaced  in  the  range  [0.0.3]ms  Relative 
to  panel  a)  the  elapsed  time  in  seconds  is  [  18,b)],  [37,c)],  [44,d)],  [55,e)],  [80,f)]. 
Note  the  expanded  y— scale  beginning  with  panel  c). 


forces  laterally  squeeze  the  interior  flow  enhancing  the  m— perturbation.  This  fur¬ 
ther  amplifies  the  initial  vertical  vorticity,  reinforces  the  vortex  force  Vjy  and  the 
instability  grows.  The  vertical  gradient  in  Stokes  drift,  acting  over  multiple  wave 
periods,  then  tilts  the  vertical  vorticity  downwind  creating  streamwise  vorticity 
cov  i  in  a  rapid-distortion  process  as  discussed  by  Leibovich  [3]  and  Teixeira  & 
Belcher  [9]. 

We  interrogated  the  LES  solutions  reported  in  [10]  looking  for  evidence  for  the 
CL2  mechanics  discussed  in  figure  7.  Classic  streamwise  oriented  Langmuir  cir¬ 
culations  are  found  in  the  presence  of  breaking  waves  similar  to  those  reported  by 
other  investigators  [13;  15]  However,  we  also  notice  that  if  breaking  occurs  at  large 
scales  and  generates  significant  amounts  of  vorticity,  then  the  convergence  result¬ 
ing  from  powerful  spanwise  vortex  forces  ±V/y  leads  to  vigorous  downwelling 
“jets”.  The  lifecycle  of  a  breaker  induced  “jet”  is  shown  in  figure  8.  The  lateral 
convergence  of  vertical  vorticity  generated  by  a  20m  breaker  under  the  action  of 
the  vortex  force  is  clearly  evident  in  this  example.  We  found  that  these  coherent 
structures,  which  persist  well  below  the  water  surface,  co-exist  with  traditional 
Langmuir  cells.  Overall  our  LES  results  suggest  that  breaker  vorticity  couples 
with  the  CL2  instability  mechanism  and  unique  coherent  structures  emerge  from 
the  interactions  between  waves  and  turbulence. 
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The  wind-driven  stably  stratified  mid-latitude  oceanic  surface  turbulent  boundary 
layer  is  computationally  simulated  in  the  presence  of  a  specified  surface  gravity-wave 
field.  The  gravity  waves  have  broad  wavenumber  and  frequency  spectra  typical  of 
measured  conditions  in  near-equilibrium  with  the  mean  wind  speed.  The  simulation 
model  is  based  on  (i)  an  asymptotic  theory  for  the  conservative  dynamical  effects 
of  waves  on  the  wave-averaged  boundary-layer  currents  and  (ii)  a  boundary-layer 
forcing  by  a  stochastic  representation  of  the  impulses  and  energy  fluxes  in  a  field  of 
breaking  waves.  The  wave  influences  are  shown  to  be  profound  on  both  the  mean 
current  profile  and  turbulent  statistics  compared  to  a  simulation  without  these  wave 
influences  and  forced  by  an  equivalent  mean  surface  stress.  As  expected  from  previous 
studies  with  partial  combinations  of  these  wave  influences,  Langmuir  circulations  due 
to  the  wave-averaged  vortex  force  make  vertical  eddy  fluxes  of  momentum  and 
material  concentration  much  more  efficient  and  non-local  (i.e.  with  negative  eddy 
viscosity  near  the  surface),  and  they  combine  with  the  breakers  to  increase  the 
turbulent  energy  and  dissipation  rate.  They  also  combine  in  an  unexpected  positive 
feedback  in  which  breaker-generated  vorticity  seeds  the  creation  of  a  new  Langmuir 
circulation  and  instigates  a  deep  strong  intermittent  downwelling  jet  that  penetrates 
through  the  boundary  layer  and  increases  the  material  entrainment  rate  at  the  base 
of  the  layer.  These  wave  effects  on  the  boundary  layer  are  greater  for  smaller  wave 
ages  and  higher  mean  wind  speeds. 


1.  Introduction 

The  air-sea  interface  and  more  broadly  the  ocean  mixed  layer  (or  ocean  boundary 
layer,  OBL)  plays  an  important  role  in  geophysical  flows.  It  modulates  air-sea  fluxes 
in  the  global  ocean  circulation  (Large,  McWilliams  &  Doney  1995;  McWilliams 
1996),  and  in  the  most  extreme  environment  controls  tropical  cyclone  intensity 
(Emanuel  2004).  The  ocean  surface  layers  support  air-sea  fluxes,  surface  gravity  waves, 
boundary-layer  turbulence,  Ekman  currents  and  significant  air  (gas)  entrainment 
when  the  surface  is  disrupted  by  intermittent  breaking  waves.  Despite  the  differences 
in  length  scales  0(1  mm-100  m)  and  time  scales  0(1  ms-1  h),  the  physics  of  these 
phenomena  are  closely  linked.  Progress  by  direct  computational  methods  has  been 
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impeded  by  the  structural  complexity  of  upper-ocean  flows,  disparity  in  magnitude 
and  time  scales  among  the  flow  components,  and  by  high  Reynolds  numbers  Re.  The 
current  generation  of  bulk  (ensemble  average)  computational  models  of  the  OBL 
used  in  large-scale  prediction  codes  does  not  account  for  the  dynamic  nature  of  the 
surface  wave  field.  These  1-D  (vertical  column)  parameterizations  of  the  OBL  (e.g. 
Mellor  &  Yamada  1982;  Large  et  al.  1995;  Price,  Weller  &  Pinkel  1986)  use  vertical 
mixing  rules  partially  based  on  past  experience  with  the  atmospheric  boundary  layer 
and  laboratory  experiments.  This  perception  of  the  OBL  has  been,  in  part,  a  practical 
necessity  as  the  hostile  measuring  environment  does  not  permit  easy  probing  of  the 
important  interactions  that  determine  the  turbulent  fluxes  of  momentum,  scalars  and 
energy  across  the  temporally  and  spatially  evolving  air-sea  interface.  While  many 
aspects  of  surface  waves  are  well  known,  a  clear  understanding  of  their  relationship 
to  fluxes,  turbulence  and  currents  has  long  been  elusive.  However,  the  impact  of  waves 
on  OBL  dynamics  is  acknowledged  in  more  recent  bulk  models  (e.g.  Craig  &  Banner 
1994;  McWilliams  &  Sullivan  2000;  Smyth  et  al.  2002;  Kantha  &  Clayson  2004; 
Lewis  &  Belcher  2004;  Melsom  &  STEtra  2004;  Ardhuin  &  Jenkins  2006;  Rascale, 
Ardhuin  &  Terray  2006),  with  improvements  motivated  by  better  measurements  and 
results  from  turbulence-resolving  simulations  with  wave  effects  as  described  here. 

The  technique  of  large-eddy  simulation  (LES)  calculates  the  currents  and  turbulence 
in  response  to  the  fluxes  using  an  embedded  model  for  subgrid-scale  turbulent  mixing 
and  dissipation.  Relatively  recently,  the  wave  influences  associated  with  Stokes  drift  uSt 
(i.e.  vortex  force  and  Lagrangian  material  transport)  have  been  incorporated  in  LES 
(Skyllingstad  &  Denbo  1995;  McWilliams,  Sullivan  &  Moeng  1997;  Polton,  Lewis  & 
Belcher  2004;  Li,  Garrett  &  Skyllingstad  2005).  These  influences  represent  conservative 
wave-current  interactions  averaged  over  the  more  rapid  wave  fluctuations  (Craik  & 
Leibovich  1976;  McWilliams  &  Restrepo  1999;  McWilliams,  Restreppo  &  Lane  2004). 
Their  effect  is  to  induce  coherent  Langmuir  circulations  in  the  OBL  turbulence  that 
importantly  modify  the  resulting  mixing,  dissipation  and  Ekman  current  profile.  In 
shallow  water,  depth-filling  Langmuir  cells  also  play  important  roles  in  sediment 
transport  (Gargett,  Tejada-Martinez  &  Grosch  2004).  Thorpe  (2004,  2005)  provides 
recent  reviews  of  Langmuir  circulations  and  their  role  in  upper  ocean  dynamics.  We 
note  that  low-Reynolds  number  direct  numerical  simulation  (DNS)  and  LES  have 
previously  been  used  to  examine  turbulent  water  flows,  namely  free-surface  flows  and 
stress-driven  interfaces,  where  the  wave  amplitudes  and/or  the  wave  phase  speeds  are 
small  compared  to  those  considered  here  (e.g.  Enstad,  Nagaosa  &  Alendal  2006;  Tsai 
1998;  Tsai  &  Yue  1996). 

Intermittent  wave  breaking  is  also  important  in  this  regime,  but  it  has  rarely  been 
included  in  oceanic  LES.  Craig  &  Banner  (1994)  represented  breaking  as  a  surface 
flux  of  turbulent  kinetic  energy  which  enhances  mixing  and  dissipation  rates  near 
the  surface  in  a  Reynolds-averaged  boundary-layer  model,  and  Noh,  Min  &  Raasch 
(2004)  verified  these  effects  in  LES.  Furthermore,  breaking  waves  transfer  momentum 
from  the  wind-generated  wave  field  to  the  oceanic  turbulence  and  currents;  this  is  the 
dominant  mechanism  for  the  net  momentum  exchange  between  winds  and  currents 
usually  associated  with  the  wave-averaged  surface  wind  stress  (Donelan  1998;  Banner 
&  Peirson  1998).  The  conventional  LES  practice  is  to  specify  a  mean  surface  stress  r. 
However,  this  does  not  adequately  represent  the  intermittent  injection  of  momentum 
and  energy  from  the  wave  field  by  breaking.  An  alternative  LES  representation  of 
breaking  that  captures  both  energy  and  momentum  injection  (Sullivan,  McWilliams 
&  Melville  2004)  is  a  random  field  of  near-surface  impulses  whose  attributes  are 
designed  to  match  field  and  laboratory  observations  of  breaking  waves  (e.g.  Rapp 
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&  Melville  1990;  Melville  &  Matusov  2002;  Melville,  Veron  &  White  2002)  and  the 
mean  fluxes  of  momentum  and  energy. 

This  paper  reports  LES  for  the  OBL  with  both  types  of  wave  effects  predominantly 
for  approximate  equilibrium,  or  ‘fully-developed’,  conditions  among  the  wind,  waves 
and  currents.  However,  we  also  show  that  the  qualitative  effects  of  the  wave  processes 
are  retained  for  smaller  wave  ages,  or  fetch-limited  conditions.  For  a  given  mean 
wind  U„  (e.g.  the  wind  at  the  height  z  =  10  m)  and  sufficiently  large  fetch,  the  wave 
field  comes  into  a  fully-developed  wavenumber-spectral  distribution  for,  say,  the  sea- 
surface  elevation,  with  generation  by  wind  form  stress  balancing  loss  by  breaking 
integrated  over  the  spectrum.  The  elevation  spectrum  implies  a  mean  Stokes  drift 
profile  iis' (c),  and  the  breaker  spectrum  implies  a  probability  distribution  function  for 
the  random  impulses.  The  impulse  distribution  is  normalized  to  have  the  equivalent 
integrated  current  acceleration  expected  from  the  empirical  bulk  formula  relating  the 
mean  wind  and  wind  stress  r,  as  well  as  an  energy  input  equivalent  to  the  empirical 
kinetic  energy  loss  rate  from  the  waves.  Thus,  the  impulses  can  wholly  replace  r  in 
the  conventional  LES  formulation  and  be  a  more  faithful  representation  of  the  true 
mechanisms  of  air-sea  momentum  and  energy  fluxes  through  the  wave  field. 

Our  results  exhibit  a  subtle  interplay  between  the  vortex  force  and  breaking. 
Langmuir  circulations  intensify  the  surface  vorticity  and  enhance  turbulent  transport 
throughout  the  boundary  layer  into  the  underlying  entrainment  zone.  Breakers 
enhance  the  near-surface  mixing,  destroy  the  Monin-Obhukov  similarity  structure, 
provide  seed  vorticity  for  wave-current  (CL2,  Leibovich  1983)  instabilities,  and  limit 
the  spatial  extent  of  the  Langmuir  circulations.  Under  certain  circumstances  breaker 
responses  and  Langmuir  circulations  combine  to  make  strong  intermittent  downward 
jets  that  further  enhance  the  vertical  transport  efficiency. 

Section  2  presents  the  representation  of  wind  and  waves  in  near-equilibrium 
conditions.  Section  3  presents  the  stochastic  breaker  model.  Section  4  is  the  LES 
model  formulation.  Section  5  describes  the  experimental  design.  Section  6  is  an 
analysis  of  a  primary  experiment  with  Ua  =  15  m  s-1.  Section  7  analyses  model 
sensitivities,  especially  to  wind  speeds  up  to  a  value  of  Ua  =  30  m  s-1,  and  wave  age. 
Finally,  §  8  is  a  summary  and  discussion. 


2.  Winds  and  waves  approaching  equilibrium 

For  modeling  wave  effects  in  the  ocean  mixed  layer,  we  must  specify  both  the  wave 
spectrum  to  determine  the  profile  of  Stokes  drift  and  the  distribution  (spectrum)  of 
breaking  waves  for  transfer  of  momentum  and  kinetic  energy  to  the  underlying  water 
column.  In  both  cases,  our  primary  interest  is  the  situation  approaching  wind-wave 
equilibrium  where  the  wave  spectrum  is  well  developed  and  the  atmospheric  inputs 
of  momentum  and  energy  flux  are  perhaps  best  characterized  by  observations.  Full 
development  refers  to  an  asymptotic  state  in  which  the  energy  density,  the  peak 
frequency  and  the  shape  of  the  wave  spectrum  are  changing  very  slowly  with  time 
following  synoptic  changes  in  the  wind.  We  neglect  the  portion  of  wind  input  which 
goes  into  bulk  wave  growth  in  the  action  balance  equation  (Komen  et  al.  1994,  p.  47) 
and  always  adopt  the  leading-order  balance  of  wind  input  equal  to  dissipation  by 
breaking.  Phillips  (2002)  has  developed  a  theoretical  model  for  the  unsteady  influence 
of  waves  on  Langmuir  circulations. 

Apart  from  overall  equilibrium  between  wind  and  waves,  or  full-development, 
an  equilibrium  subrange  in  a  developing  wind-wave  spectrum  may  also  exist.  This 
subrange  corresponds  to  the  region  in  which  the  sum  of  the  wind  input,  nonlinear 
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wave-wave  interactions  and  the  dissipation,  mainly  due  to  breaking,  are  in  balance 
(Phillips  1985). 

2.1.  Wave  spectra  and  Stokes  drift 

The  vertical  profile  of  Stokes  drift  us'(z)  is  a  primary  input  to  OBL  models 
that  incorporate  wave-averaged  effects  (e.g.  McWilliams  et  al.  1997;  Kantha  & 
Clayson  2004).  Past  LES  typically  estimate  us,(z)  from  monochromatic  wave  fields 
(Skyllingstad  &  Denbo  1995;  McWilliams  et  al.  1997;  Noh  et  al.  2004)  chosen  to 
approximate  a  dominant  component  in  the  wave  spectrum  (the  work  described 
by  Polton  et  al.  (2004)  is  an  exception).  In  the  case  of  wind-wave  equilibrium, 
measurements  of  the  wave  spectrum  are  more  complete  and  thus  we  can  be  more 
precise  with  the  specification  of  the  Stokes  profile.  Recently,  Alves,  Banner  &  Young 
(2003)  revisited  the  classical  Pierson-Moskowitz  wave-height  spectrum  (Pierson  & 
Moskowitz  1964)  for  fully-developed  seas,  reanalysing  the  original  database  to  refine 
the  spectral  shape  and  the  integral  parameters  that  quantify  the  spectrum.  A  database 
of  29  events  was  interrogated  to  establish  the  spectral  shape  and  fitting  constants. 
Adhering  to  the  spectral  shape  proposed  by  Donelan,  Hamilton  &  Hui  (1985),  Alves 
et  al.  (2003)  favor  the  empirical  spectral  form, 


F(f) 


awg2 

(2  7X)4  /4  /„exp 


(2.1) 


where  /  =  a/2n  is  the  frequency,  er  is  the  radial  frequency,  aw  is  the  ‘Phillips  constant’, 
g  is  the  gravitational  acceleration  and  fp  is  the  peak  in  the  spectrum  related  to  a 
reference  atmospheric  wind  Ua(z  =  10  m)  by 

v  =  fPUa/g.  (2.2) 

Empirical  constants  that  appear  in  the  above  expressions  are  v  «  0.123  and  aw  « 
6.15  x  10~3.  The  observations  of  significant  wave  height,  Hs  =  4  (iff12,  where  (if)  is 
the  variance  of  the  surface  displacement  due  to  the  waves,  closely  follow  the  curve  fit 

Hs  =  0.24  Ul/g.  (2.3) 

The  wave  energy  density  E  =  p0g(ri2),  where  pa  is  the  density  of  the  water.  The  phase 
speed  at  the  peak  of  the  equilibrium  wave  height  spectrum  is  cPiE  =  Ua/2tiv.  Thus 
the  phase  speed  at  the  peak  is  moving  about  20%  faster  than  the  reference  wind  in 
fully-developed  conditions.  The  linear  dispersion  relationship  c2  =  g/k  along  with 
/  =  g/2nc  is  used  to  move  between  frequency  /,  wavenumber  k  and  phase  speed 
c  spaces.  The  observations  that  are  used  to  determine  the  wave  spectral  shape  and 
constants  extend  out  to  f/fp  »  4. 

McWilliams  &  Restrepo  (1999,  equation  62),  using  the  procedure  outlined  by 
Kenyon  (1969),  show  how  to  incorporate  a  full  wave  spectrum  in  the  estimate  of 
uSt.  For  our  LES  (§4),  we  follow  this  approach  and  employ  numerical  quadrature  to 
evaluate 

us'(z)  =  —  /  F(er)crexp  - — 4  dcr .  (2.4) 

g  Jo  L  g  . 

The  main  difference  between  the  Alves-Donelan  wave  spectrum,  given  by  (2.1), 
and  the  Pierson-Moskowitz  formula  is  that  the  former  decays  at  a  slower  rate  at 
high  frequencies,  f~4  compared  to  /~5.  This  slope  change  has  no  influence  on  uSt 
for  z  <  —2  m,  but  for  z  =  [—2,  0]  m  the  Alves-Donelan  spectrum  generates  a 
sharper  near-surface  gradient  in  the  Stokes  profile  compared  to  that  obtained  from 
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the  Pierson-Moskowitz  wave  spectrum.  A  high-frequency  tail  with  slope  f~5  can  be 
appended  to  the  spectral  shape  given  by  (2.1)  to  extend  its  range  of  applicability  (e.g. 
Banner  1990;  Terray  et  al.  1996). 

2.2.  Bulk  atmospheric  inputs  of  momentum  and  energy 
At  wind  speeds  above  (6-8)  m  s-1,  the  wave  field  is  the  main  transmitter  of 
atmospheric  momentum  and  energy  passing  nearly  100%  of  these  fluxes  locally 
to  the  underlying  currents  (Terray  et  al.  1996;  Donelan  1998,  2001)  primarily  through 
the  action  of  wave  breaking.  Thus  the  global  amount  of  momentum  and  energy  that 
is  available  to  generate  currents  from  breaking  is  known  based  on  the  atmospheric 
inputs.  The  bulk  aerodynamic  method  is  most  often  used  to  determine  the  momentum 
flux  in  the  atmospheric  surface  layer.  There  is  debate  about  the  empirical  constants 
appearing  in  the  formula 

Xa=  PaCd\Ua\Ua,  (2.5) 

but  this  functional  form  is  generally  well  accepted.  Here  pa  is  the  air  density,  Ua  is 
the  mean  wind  at  the  reference  height  z  =  10  m,  and  Cd  is  a  wind  speed  dependent 
drag  coefficient.  The  atmospheric  friction  velocity  is  then  u»a  =  y/\xa\/pa-  For  our 
purposes  we  adopt  the  drag  parameterization  of  Liu,  Katsaros  &  Businger  (1979), 


Cd  = 


1.3  x  10 


,-3 


(0.79  +  0.0509  x  Ua)  x  10 


i-3 


Ua  ^  10  m  s  1 
Ua^  10  m  s-1 


(2.6) 


noting  other  prescriptions  are  available  (e.g.  Large  &  Pond  1982;  Donelan  1998; 
Fairall  et  al.  2003)  that  include  the  effects  of  finite  fetch  and  the  state  of  wave 
development  (or  wave  age).  Observations  indicate  that  the  linear  variation  of  Cd 
with  wind  speed  continues  up  to  about  25-30  m  s-1.  For  even  higher  winds  there 
are  few  field  measurements  (Powell,  Vickery  &  Reinhold  2003),  but  both  field  and 
laboratory  data  show  that  the  drag  coefficient  saturates,  or  reaches  a  weak  maximum, 
around  Cd  *  2.3  x  10~3  at  wind  speeds  around  30-35  m  s_1  (Donelan  et  al.  2004). 
The  causes  of  this  regime  change  in  Cd  near  the  threshold  for  hurricane  winds  are  a 
topic  of  current  research  since  the  maximum  intensity  of  hurricanes  depends  on  the 
ratio  of  the  coefficients  of  momentum  to  enthalpy  transfer  across  the  air-sea  interface 
(Emanuel  2004). 

Breaking  waves  supply  momentum  for  current  generation  but  also  serve  as  a 
significant  source  of  kinetic  energy  for  the  water  column.  Our  wave-OBL  coupling 
accounts  for  the  energy  flux  flowing  from  winds  to  waves  to  currents.  There  are  no 
direct  measurements  of  energy  flux  from  the  atmosphere  to  the  wave  field  and 
thus  an  alternative  method  is  required  to  determine  this  variable.  Terray  et  al.  (1996) 
propose  the  parameterization 


£«  =  K\ c 


where  c  =  u,a  g,(cp/u,a) , 


(2.7) 


based  on  a  wave-growth  model  developed  by  Donelan  &  Pierson  (1987).  In  this 
expression  c  is  an  ‘effective  phase  speed’,  a  velocity  scale  characteristic  of  the  wave 
field,  and  it  reflects  contributions  over  the  entire  wave  spectrum.  Based  on  six  different 
observational  datasets,  Terray  et  al.  (1996,  figure  8)  find  that  c/u,a  (or  gt)  depends 
critically  on  the  state  of  wave  development.  We  adopt  wave  age,  defined  as  the  ratio 
of  the  phase  speed  at  the  peak  of  the  wave  spectrum  to  the  atmospheric  friction 
velocity  cp/u,a,  as  a  bulk  measure  of  wave  development,  g,  is  in  the  range  3-4  for 
waves  approaching  full  development  ( cp/u,a  «  30),  attains  a  maximum  of  about  7-9 
for  developing  seas  ( cp/u,a  =  15-20),  and  falls  to  small  values  approximately  1-4  for 
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very  young  waves  (cp/u»a  <  10).  The  explicit  conversion  between  wave  age  cp/u,a 
and  Terray  parameter  g,  used  in  the  present  work  is  given  in  table  1  in  §  5.  Terray 
et  al.  (1996)  find  that  dissipation  measurements  taken  beneath  breaking  waves,  when 
scaled  by  the  wind  input  £„  and  significant  wave  height  Hs,  collapse  their  data  over 
a  range  of  growing  waves  4.3  <  cp/u,a  <  1  A.  Further  support  for  the  energy  flux 
parameterization  (2.7)  is  provided  by  Drennan  et  al.  (1996)  as  they  show  this  form 
also  appears  to  apply  for  larger  wave  ages  closer  to  wind-wave  equilibrium.  Phillips 
(1985)  (see  also  Melville  1994)  has  developed  an  equilibrium  model  for  atmospheric 
energy  flux  based  on  breaking  statistics  that  is  functionally  different  from  the  proposal 
by  Terray  et  al.  (1996),  but  its  predictions  are  in  general  agreement  with  (2.7).  In 
the  absence  of  additional  information,  we  adopt  (2.7)  as  our  prescription  of  the 
atmospheric  energy  flux  and  consider  the  impacts  of  wave  age  in  a  limited  range 
approaching  equilibrium  cp/u,a  =  (19-30). 

As  a  consistency  check  we  note  that  in  the  absence  of  other  external  energy  sources 
the  vertically  integrated  dissipation  rate  in  the  water  is  approximately  equal  to  the 
average  energy  flux  £fl  from  winds  to  waves: 

£a/Po  »  f  edz,  (2.8) 

J  —00 

where  e  is  the  local  dissipation  rate  per  unit  of  water  mass  and  z  extends  over 
the  depth  of  the  water  column.  This  integral  balance  is  in  accord  with  the  basic 
assumption  that  breaking  is  the  dominant  energy  source  (Terray  et  al.  1996).  In  our 
simulations  we  use  (2.8)  as  an  indicator  of  flow  stationarity  and  global  conservation  of 
energy  between  the  atmospheric  inputs  and  dissipation  in  the  water  column.  However, 
the  LES  also  includes  current  shears,  stratification  and  Stokes  drift  effects  that  further 
increase  the  dissipation  in  the  water  in  addition  to  the  contributions  from  breaking 
waves. 

At  this  point,  notice  that  the  specification  of  the  wave  fields,  Stokes  drift  and  fluxes 
of  momentum  and  energy  can  all  be  expressed  in  terms  of  the  reference  atmospheric 
wind.  Hence  the  primary  atmospheric  input  to  our  LES  model  of  the  OBL  is  Ua. 
Developing  seas  are  parameterized  by  wave  age  as  well  as  wind  speed. 

3.  Stochastic  breakers 

The  other  physical  process  we  wish  to  capture  in  our  modeling  is  the  intermittent 
and  dynamical  influences  of  breaking  waves.  Breaking  encompasses  a  wide  range 
of  scales  and  is  sufficiently  complex  that  we  must  idealize  this  process,  retaining 
its  essential  ingredients.  First,  we  assume  that  the  onset  of  breaking-which  is  likely 
linked  to  the  local  winds,  wave-wave  and  wave-current  interactions  (Melville  1996)— is 
a  random  process  uniformly  distributed  in  space  and  time.  This  is  consistent  with  the 
visual  appearance  of  the  broken  sea  surface  under  high  wind  conditions.  However, 
it  is  also  well  known  that  the  incidence  of  large-scale  breaking  is  related  to  the 
group  structure  of  the  surface  wave  field,  in  an,  as  yet,  not  fully  understood  way 
(Donelan,  Longuet-Higgins  &  Turner  1972).  In  the  present  formulation  the  averaging 
is  over  time  scales  much  larger  than  the  characteristic  time  of  the  wave  groups,  so 
any  departure  from  an  assumption  of  random  incidence  of  breaking  is  not  justified. 
Secondly,  we  represent  an  individual  breaker  as  a  compact  function  in  space  and 
time  that  captures  the  (average)  bulk  impulse  from  breaking  (Sullivan  et  al.  2004). 
This  avoids  the  considerable  complexity  of  a  full  air-water  simulation  that  cannot 
simultaneously  span  the  range  of  scales  from  bubbles  to  mixed-layer  turbulence.  Also, 
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the  imposed  breaking  impulses  are  assumed  to  be  non-interacting,  which  allows  a 
field  of  breakers  to  be  built  from  a  linear  superposition  of  discrete  events. 

3.1.  Momentum  transfer  from  waves  to  currents 
The  wave  field  is  an  intermediate  reservoir  that  releases  momentum  to  the  underlying 
water  in  localized  intermittent  impulses,  and  the  primary  path  to  current  generation 
is  through  wave  breaking.  The  momentum  passed  to  the  OBL  is  then  a  summation  of 
randomly  occurring  breaking  events.  In  our  modeling  the  acceleration  from  a  single 
event  has  only  horizontal  components, 

An  =  t,  cm)(jccos  0m  +  y  sin  0m) ,  (3.1) 


where  0m  is  the  breaker’s  angle  of  orientation  in  the  LES  coordinate  system  and  sdm 
is  our  parameterization  of  the  local  space-time  structure  of  acceleration  amplitude 
for  a  single  event  in  terms  of  its  phase  speed,  cm,  and  initiating  location  and  time, 
(jc "ft'")  (§3.2).  The  breaker’s  phase  speed  vector  is  d"  =  c"'( it  cos©'”  +  ysin@m).  a 
single  event  provides  momentum 

fT(cm)  r  V(cm ) 

Jfn(cm,  ©'”)  =  Po  \  /  Am( JC,  t,  cm,  0m)  dx  At  (3.2) 

Jo  Jo 

with  magnitude 

rn<t")  rV(cm) 

\Jtm\  =Po  /  t,cm)  dx  dr .  (3.3) 

Jo  Jo 


The  integration  limits  reflect  the  compact  duration  of  breaking  over  space  and  time 
relative  to  the  initiating  coordinates,  i.e.  the  impulse  is  non-zero  over  the  period,  T(c), 
and  volume,  V  (c),  of  the  breaking  wave. 

The  total  momentum  supplied  to  the  currents  is  the  integrated  effect  of  many 
randomly  occurring  breaking  events  of  varying  size,  speed  and  orientation, 


Mb 


Jt(c,  0)c0*(c,  0)AcA0  . 


(3.4) 


In  this  expression,  N  is  the  total  number  of  breakers  and  3P(c,  0)  is  the  probability 
density  function  (PDF)  of  breaker  speeds  c  with  orientation  0  and  is  normalized 
according  to  the  definitions, 


P(c)  = 


-f 


c^>(c,  0)d0,  where 


P(c)  dc  =  1 . 


(3.5) 


N  and  0*{c,  0)  are  functions  of  the  environmental  conditions,  in  particular  wind 
speed  and  wave  age  (e.g.  Wu  1988;  Gemmrich  &  Farmer  1999;  Melville  &  Matusov 
2002). 

Momentum  conservation  between  atmosphere  and  ocean  limits  the  total  amount 
of  momentum  that  can  be  passed  through  the  wave  field.  The  long-time,  large-area 
total  momentum  conservation  rule  expressed  in  terms  of  the  horizontal  momentum 
flux,  ra,  from  the  atmosphere  over  an  area  of  surface  water,  and  time  period,  2T  p, 
is 


TaPPp  =  Mb,  (3.6) 

where  Mh  is  the  momentum  transferred  by  a  field  of  breakers.  Substitution  of  (3.4) 
into  (3.6)  connects  the  bulk  atmospheric  momentum  to  the  probability  of  breaking 
and  also  defines  a  critical  breaking  parameter  in  the  FES.  We  define  the  average  rate 
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of  breaker  generation  as 


N 


N 


(3.7) 


where  N  is  the  number  of  breakers  created  per  unit  of  water  surface  area  per  unit 
of  time  (units  of  s_1  nr2).  In  other  words,  a  sufficient  number  of  breakers  N  at  the 
proper  scale  must  be  created  at  a  rate  N  to  match  the  long-time  large-area  average 
of  momentum  (and  energy)  transfer  in  air  and  water  given  by  (3.6).  It  is  expected 
that  N  =  N(Ua ,  cp/u,a). 

Breaking  releases  momentum  to  generate  currents  but  also  energizes  the  water 
column,  which  is  vividly  apparent  from  images  of  breaking  in  the  laboratory  and  in 
the  field.  Thus  our  modeling  also  needs  to  account  for  energy  transfer  from  winds 
to  waves  to  currents.  From  consideration  of  the  reduction  in  momentum  and  energy 
densities  of  the  wave  field  due  to  breaking  (Phillips  1977),  the  total  (kinetic  and 
potential)  energy  released  by  an  event  is  (Phillips  1985).  This  is  consistent 

with  the  measurements  of  Terray  et  al.  (1996),  is  similar  to  that  used  in  second- 
order  closure  parameterizations  (Craig  &  Banner  1994)  and  is  consistent  with  our 
parameterization  for  breaker  momentumf.  By  analogy  with  (3.4)  and  (3.6)  the  total 
energy  transferred  by  a  field  of  breakers  is  then 


Eh  = 


c-J((c,  0)c3P(c,  0)dcd©  . 


(3.8) 


Matching  the  atmospheric  energy  flux  from  winds  to  waves  requires 

Pa<f2Tp  =  Eb.  (3.9) 

3.2.  Momentum  impulse  from  a  single  breaker 

The  breaker  model  used  here  (Sullivan  et  al.  2004)  is  based  on  laboratory 
measurements  by  Rapp  &  Melville  (1990),  Loewen  &  Melville  (1990),  Melville  et  al. 
(2002)  and  video  imagery  of  breaking  in  the  field  (Melville  &  Matusov  2002).  The 
main  results  of  the  laboratory  experiments  are  that:  the  duration  of  active  breaking 
is  0(T );  the  maximum  fluid  velocity  during  breaking  is  O(c);  and  the  depth  of 
penetration  of  the  broken  fluid  is  <9 (2a),  where  T,  c  and  a,  are  the  characteristic 
period,  linear  phase  speed  and  amplitude,  respectively,  of  the  breaking  wave.  The 
laboratory  experiments  also  show  that  the  post-breaking  flow  variables  scale  with 
the  pre-breaking  variables  even  at  large  times.  We  refer  to  this  set  of  dimensional 
estimates  as  Froude  scaling  of  breaking,  with  length  and  time  scales  related  through 
the  linear  dispersion  relationship.  It  is  the  basis  of  our  formulation  of  the  breaking 
model. 

It  may  seem  curious  that  a  strongly  nonlinear  process  such  as  breaking  can  be 
scaled  using  a  linear  dispersion  relationship.  However,  there  are  some  rather  strong 
constraints  on  the  theoretical  maximum  and  empirical  maximum  phase  velocities  for 
surface  gravity  waves,  and  it  is  the  phase  velocity  that  relates  the  length  and  time 
scales.  The  maximum  phase  speed  for  nonlinear  surface  gravity  waves  is  9%  larger 
than  the  linear  phase  speed  and  occurs  at  a  wave  slope  of  0.436,  slightly  less  than 
the  slope  of  0.443  for  the  Stokes  limiting  wave  (Longuet-Higgins  1975).  However, 
in  practice  fast  instabilities  and  breaking  tend  to  limit  the  maximum  practical  wave 


f  At  low  winds  breaking  is  less  obviously  the  dominant  air-sea  transfer  process.  Our  breaker 
representation  is  thus  relevant  to  the  wind  regime  with  Ua  greater  than  about  5  m  s_1  (Makin, 
Kudryavtsev  &  Mastenbroek  1995;  Banner  &  Peirson  1998). 
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slope  to  0.3,  for  which  the  corresponding  steady  wave  phase  velocity  is  just  5% 
larger  than  the  linear  phase  speed.  Thus  for  breaking  kinematics,  the  use  of  the  linear 
phase  speed  may  involve  errors  of  around  5%.  For  the  dynamics,  where  the  breaking 
momentum  and  energy  fluxes  scale  as  the  fourth  and  fifth  powers  of  the  phase  speed 
(Phillips  1985),  the  errors  for  individual  breaking  events  may  be  up  to  22%  and  28%, 
respectively.  However,  the  number  of  breakers  is  constrained  to  satisfy  the  bulk  fluxes, 
so  in  this  model,  errors  will  flow  through  to  the  breaking  statistics.  Given  the  order 
of  magnitude  scatter  in  observations  of  breaking  statistics,  any  errors  in  the  scaling 
are  within  the  scatter  of  the  current  observations. 

The  formula  for  an  individual  breaker  impulse  is 


/•GO 

S=  s/d  t,  (3.10) 

Jo 


where 

s/  =  kb^~  FT(a) X{p)<Sf(S)Z(y),  (3.11) 

2k 

with  the  (FT,  9C ,  ,  2£)  space-time  shape  functions  that  are  determined  from  DNS  of 

a  single  isolated  breaker  designed  to  replicate  the  laboratory  results  of  Melville  et  al. 
(2002).  All  breakers  are  assumed  to  be  self-similar  and  separable  in  dimensionless 
time  and  space  coordinates, 


a  = 


P 


c(t-t0 )  ’ 


=  2(y  -  y0) 
A 


z 

xc{t  -  to)  ’ 


(3.12) 


with  (t0,  x0,  y„,  z0  =  0)  the  onset  time  and  position  of  the  chosen  breaker;  (c,  A,  T) 
are  its  phase  speed,  wavelength  and  period.  The  wave  characteristics  (c,A,T)  are 
related  to  each  other  through  the  deep-water  linear  dispersion  relation  c2  =  gAf 2n 
with  T  =  A/c.  Notice  that  in  the  normalized  vertical  coordinate  (3.12)  the  constant 
0  <  x  <  U  which  is  just  the  aspect  ratio  of  the  depth  to  length  of  the  breaker, 
controls  the  depth  penetration  of  the  breaker  forcing;  x  =  0-2  matches  the  laboratory 
measurements  of  Melville  et  al.  (2002).  The  specific  functional  forms  for  the  shape 
functions  is  the  set  of  equations  (3.3)  given  in  Sullivan  et  al.  (2004).  The  momentum 
supplied  by  our  breaking  model  (3.11)  varies  with  wave  phase  speed  c  because  of  the 
explicit  dependences  in  (3.12).  If  we  substitute  the  model  (3.11)  for  s/  into  (3.3)  for 
breaker  momentum  and  convert  the  time  and  space  integrations  dx  dr  — *  dp  d<5  d y  da 
we  find  the  breaker  momentum  increases  rapidly  with  phase  speed,  i.e.  Jt(c)  ~  c1 . 
Similarly  the  breaker  energy  grows  as  c8.  These  power-law  dependences  impact  the 
determination  of  the  breaker  PDF. 


3.3.  Breaker  PDF 

Bulk  conservation  of  momentum  and  energy  fluxes  between  the  atmosphere  and  ocean 
constrains  the  properties  of  the  breaker  field.  Given  Ea  Af:Tp  =  Eh  and  xa  Af.Tp  =  Mh 
and  atmospheric  inputs  the  unknowns  are  then  the  distribution  of  breakers  across  c 
(here  the  PDF  of  c)  and  the  breaker  generation  rate.  N  is  first  eliminated  by  forming 
the  energy-momentum  flux  ratio, 

£a/|Tj  =  Eb/\Mh\.  (3.13) 


Substitution  of  the  parameterization  for  atmospheric  energy  flux,  given  by  (2.7),  and 
the  definitions  of  breaker  momentum  and  energy,  given  by  (3.4)  and  (3.8),  into  (3.13) 
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M*a  §t 


®)|  cZP(c,  0)  dc  d® 


c-Jt{c ,  0)c^(c,  ®)dcd®  . 


(3.14) 

Equation  (3.14)  is  a  weak  global  constraint  on  the  breaker  PDF.  The  distribution 
of  (discrete)  breakers  with  phase  speed  and  direction  is  unknown  and  additional 
information  is  needed  to  develop  a  specific  computational  form  for  P(c ,  0). 
Quantifying  the  distribution  and  magnitude  of  breakers  in  the  field  is  a  challenging 
observational  task  and  as  a  result  the  wave-breaking  dissipation  spectrum  remains 
the  least  understood  term  in  the  wave-action  equation  used  in  numerical  wave 
forecasting  (Komen  et  al.  1994;  Donelan  2001).  Recently,  Melville  &  Matusov  (2002) 
and  Gemmrich  (2005)  attempted  to  measure  breaking  statistics  using  video  imagery 
from  an  aircraft  and  a  floating  instrument  platform  (R/P  FLIP),  respectively.  These 
measurements  were  partly  motivated  by  the  theoretical  work  of  Phillips  (1985)  that 
proposes  breaking  statistics  can  be  built  from  A(c)dc,  the  average  length  of  breaking 
crests  per  unit  area  of  ocean  surface  traveling  at  velocities  in  the  range  (c,  c  +  dc). 
Digitized  video  output  supplies  data  on  the  distribution  of  A(c),  moments  of  /1(c) 
(corresponding  to  momentum  and  energy  flux  from  waves  to  currents),  directional 
distribution  of  breakers,  as  well  as  other  statistics.  The  video  measurement  technique 
relies  on  imaging  active  whitecaps:  breakers  with  visible  air  entrainment.  Breakers 
generating  few  bubbles  or  foam  are  missed  by  the  detection  algorithms  leading  to 
an  underestimate  of  small-scale  breaking.  Melville  &  Matusov  (2002)  estimate  the 
small-scale  resolution  in  their  observations  to  be  X  «  0.7  m  corresponding  to  a  phase 
speed  c  «  1  m  s_1.  An  important  result  from  this  dataset  shows  an  exponential 
decrease  of  /1(c)  with  c  and  a  variation  of  /1(c)  as  wind  speed  cubed  over  the  range 
of  wind  speeds  7.2-13.6  m  s-1. 

A  direct  connection  between  the  PDF  of  breaking  that  appears  in  (3.14)  and 
measurements  of  /1(c)  depends  on  multiple  unknown  factors,  including  the  model  for 
breaker  momentum  transfer.  Given  this  ambiguity  we  make  two  assumptions:  (I)  we 
neglect  the  directional  distribution  of  breakers  and  assume  alignment  between  the 
surface  winds  and  breakers  noting  that  this  dependence  can  easily  be  introduced  when 
more  observational  data  become  available;  (II)  we  adopt  an  exponential  functional 
form  for  the  breaker  PDF  by  analogy  with  the  observations  of  Melville  &  Matusov 
(2002)  with  the  matching  condition  (3.14)  used  to  determine  any  unknown  coefficients. 
An  important  advantage  of  this  approach  is  that  bulk  conservation  of  momentum 
and  energy  is  enforced.  The  wind-speed-dependent  breaker  PDF  used  here  is  thus 


P(c)  =  exp(— Z?2  c/u,a). 


(3.15) 


where  (£>i,h2)  are  modeling  constants.  Our  preference  for  using  a  wind  speed 
dependence  based  on  u,a  instead  of  Ua  follows  from  the  explicit  appearance  of 
friction  velocity  in  (3.14)  and  the  numerical  evaluation  of  (3.14)  discussed  near  (3.16). 
hi  is  chosen  to  define  a  unity  PDF,  i.e.  f  P(c)dc  =  1,  while  b2  is  introduced  so  as  to 
satisfy  the  expression  given  by  (3.14). 

As  noted  above,  measurements  of  the  distribution  of  breakers  and  their  dynamical 
importance  are  incomplete  at  low  values  of  c  (i.e.  at  small  scales).  The  theoretical 
model  of  Phillips  (1985,  p.  524)  predicts  that  the  total  momentum  flux  from  wave 
breaking  depends  on  the  high  wavenumbers  at  the  upper  end  of  the  equilibrium  range 
k\.  In  this  theory,  momentum  conservation  requires  the  equilibrium  range  to  terminate 
at  a  cutoff  wavenumber  k\  =  rg/ula  or  minimum  phase  speed  ci/u,a  =  1  / ^fr  with 
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0  <  r  <  1.  Phillips  (1985)  estimates  r1/2  lies  below  0.4  but  values  as  large  as 
0.7  are  possible.  A  typical  range  of  values  for  wind  speed  Ua  =  15  m  s_1  is  then 
ci  «  0.8— 1.5  m  s-1.  Belcher  &  Vassilicos  (1997)  and  Hara  &  Belcher  (2002)  offer 
improvements  to  the  model  of  Phillips  (1985)  by  including  the  sheltering  effect  of 
long  waves  on  short  waves  in  the  equilibrium  range.  The  model  accounts  for  the 
momentum  flux  from  the  wind  to  the  longer  waves  leaving  a  smaller  stress  to  force 
growth  of  the  shorter  waves.  Assuming  that  the  dissipation  spectrum  is  proportional 
to  the  wind  input,  their  model  then  predicts  that  breaking  crests  of  shorter  waves  are 
less  frequent  especially  in  high  winds  u,a  >  0.5  m  s-1.  An  analysis  of  measured  wave 
height  spectra  by  Donelan  (2001)  also  shows  the  growth  of  small  waves  is  modulated 
by  long  waves  leading  to  a  net  reduction  in  their  energy  density,  also  implying  a 
reduced  contribution  to  fluxes  from  small  scale  breaking.  (However,  the  full  problem 
of  long-wave-short-wave  interaction,  both  through  direct  hydrodynamic  interactions 
and  indirectly  through  coupling  with  the  wind,  remains  one  of  the  unsolved  problems 
in  air-sea  interaction  and  ocean  remote  sensing.) 

Based  on  these  theoretical  predictions  and  sparse  observations  we  need  to  limit 
the  breaking  at  small  scales.  However,  given  the  uncertain  nature  of  the  available 
results,  we  simply  choose  to  truncate  the  PDF  of  breaking  at  a  small  phase  speed  c\ 
according  to  Phillips  (1985)  that  depends  on  wind  stress.  Then  (3.14)  reduces  to  the 
computational  rule, 


c/u»a]P[c)Ac  =  0, 


(3.16) 


with  the  PDF  given  by  (3.15).  The  PDF  of  breaking  depends  on  our  model  for 
breaker  momentum  flux,  wave  age  contained  in  the  Terray  parameter  g,  =  g,(cp/u,a) 
and  wind  stress.  Because  of  the  exponential  decay  of  the  PDF  at  large  c  the  integral 
converges  rapidly  and  is  only  weakly  dependent  on  the  upper  limit  of  integration 
cp.  A  rapid  decay  of  the  PDF  with  increasing  c  is  required  in  order  to  satisfy  (3.16) 
since  the  breaker  momentum  J/(c )  ~  c1 .  This  makes  the  determination  of  b2  delicate. 
Based  on  a  numerical  evaluation  of  (3.16)  we  find  for  a  given  wave  age  cp/u,a  (or 
Terray  parameter  g,)  that  a  single  value  of  the  modeling  constant  b2  is  capable 
of  satisfying  (3.16)  to  a  good  approximation  across  a  wide  range  of  wind  speeds 
(5-30  m  s-1).  This  result  supports  our  choice  of  wind  speed  scaling  based  on  u,a 
in  (3.15). 

Once  the  PDF  of  breaking  is  determined  we  retrace  our  steps  and  next  find  the 
number  of  breakers  or  breaker  generation  rate  N  by  evaluating  the  momentum 
conservation  rule  (3.6);  its  right-hand  side  is  given  by  (3.4).  Specifically,  N  is  chosen 
to  satisfy 

rcp  rT(c)  rv(c) 

xa  =  Np0  /  /  P(c)  s/(x,  t,  c)  dx  dt  dc  (3.17) 

Jci  Jo  Jo 

for  each  wind  speed  and  wave  age.  Typical  PDFs  and  breaker  generation  rates  are 
shown  in  figures  1  and  2  for  varying  Ua  and  cp/u,„.  Notice  how  the  PDF  in  figure 
1  shifts  towards  larger  scales,  larger  values  of  c,  with  increasing  winds.  Also,  for 
a  fixed  wind  speed  larger  breakers  are  obtained  for  wave  ages  in  the  intermediate 
range  cp/u,  *  19—30  (developing  seas)  compared  to  those  in  the  equilibrium  range 
cp/u,a  >  30  (fully-developed  seas).  A  similar  dependence  of  the  breaking  PDF 
on  wave  age  is  also  found  in  the  observations  of  Gemmrich  (2005,  figure  3).  The 
breaker  generation  rate  decays  as  the  wind  speed  increases  or  wave  age  decreases,  a 
consequence  of  an  overall  shift  in  the  PDF  towards  large-scale  breaking.  An  important 
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Figure  1.  The  PDF  of  breaking  P(c)  =  e~b2Ciu'a  for  varying  wind  speeds  and  small  variations 
in  wave  age.  In  (a)  P(c)  is  shown  for  wind  speeds  Ua  =  (5, 10,  15, 20, 30)  m  s-1  for 
waves  approaching  equilibrium  cp/u,a  =  30,  Phillips  parameter  r  =  0.2,  and  PDF  constant 
bi  =  2.196.  In  ( b )  the  winds  are  held  constant  Ua  =  15  m  s-1  and  the  wave  age  varies, 
Cp/u*a  =  (30,  23,  19)  while  b2  =  (2.196, 1.46,  1.06). 


attribute  of  the  exponential  PDF  is  that  it  leads  to  rapid  convergence  of  the  integrals 
(3.4)  and  (3.8)  for  total  breaker  momentum  and  energy  as  the  phase  speed  c  increases. 
Computationally,  this  means  that  a  finite  number  of  draws  from  the  PDF  will  closely 
satisfy  the  momentum  and  energy  flux  balances  between  the  atmosphere  and  ocean. 
Then  the  details  of  the  large-c  tail  of  the  PDF  are  not  critical.  We  find  more  than 
98%  of  the  breaker  momentum  and  energy  is  captured  by  breakers  satisfying  the 
criterion  P{c)  >  10~6.  Other  functional  forms  for  the  breaker  PDF  do  not  satisfy  this 
criterion  (see  Sullivan,  McWilliams  &  Melville  2005). 

The  density  distributions, 

^Jl(c,0)CP(c,0)  and  ffc-Jl(c,0)cP(c,0)  '  (318) 

provide  information  about  the  partitioning  of  breaker  momentum  and  breaker  energy 
across  the  range  of  phase  speeds,  c.  These  density  distributions  (or  breaker  spectra) 
are  normalized  so  that  their  integrals  over  (c,  0)  equal  unity  and  are  the  companions 
to  the  wave-height  spectra  given  by  (2.1).  Figure  3  shows  the  phase  speeds  (cz,cE), 
where  (3.18)  attain  their  maximum  values  for  varying  wind  speed  and  wave  age.  Notice 
when  cp/u,a  ^  30  (or  g,  >  4),  the  breaker  momentum  and  energy  distributions  peak 
at  an  intermediate  phase  speed.  For  example,  with  Ua  =  15  m  s-1  and  cp/u,a  =  23 
the  maximum  contributions  to  breaker  momentum  and  energy,  weighted  by  the  PDF, 
occurs  in  the  range  c  «  3-4  m  s-1.  Then  the  lower  limit  ci  plays  a  minor  role  in 
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Figure  2.  The  rate  of  breaker  generation,  N,  for  varying  winds  and  wave  age.  For 
( U a ,  cp / u*a)  =(15  m  s_1,  30)  and  an  x-y  spatial  domain  of  300  m  x  300  m  approximately 
5  x  106  breakers  are  generated  per  hour. 


determining  the  PDF;  the  breaker  PDF  is  more  dependent  on  wave  age  and  wind 
speed  than  on  the  minimum  phase  speed  c\.  Melville  &  Matusov  (2002)  also  found 
the  moments  c4A(c)  and  c5A(c)  reach  a  maximum  in  a  similar  range  of  phase  speeds. 
They  estimate  the  wave  age  of  their  observations  to  be  in  the  range  cp/u,a  =  16—24, 
which  corresponds  to  g,  =  8— 5. 

The  stochastic  breaking  model  described  above  can  be  used  to  make  an  estimate  of 
the  fractional  area  of  water  covered  by  breakers,  i.e.  the  whitecap  coverage  iVc,  which 
is  a  quantity  obtained  from  video  imagery.  Field  observations  of  visible  whitecaps,  i.e. 
those  with  lm  report  IV  c  ~  1.75  Uf15  over  the  wind  speed  range  JJa  =  2,  20  m  s_1 

(Wu  1983).  At  higher  wind  speeds,  closer  to  hurricane  conditions,  the  available  data 
are  limited  but  the  whitecap  coverage  is  expected  to  reach  a  limiting  value.  Predictions 
from  the  rectilinear  breaking  model  described  and  the  breaker  PDF  shown  in 
figure  1  with  a  lower  bound  c  ^  1.25  m  s-1  (or  A  =  1  m)  are  in  good  agreement  with 
this  empirical  correlation,  but  saturate  at  a  wind  speed  near  Ua  ~  17  m  s-1.  This 
agreement  shows  the  Froude  scaling  used  to  build  the  discrete  breaker  model  holds 
over  a  wide  range  of  wind  speed. 

The  dependence  of  P{c)  and  N  on  wind  speed  and  wave  age,  shown  in  figures  1  and 
2,  has  important  implications  for  the  interactions  between  vortex  force  and  breaking 
as  discussed  in  §  6.  For  a  given  wind  speed,  growing  seas  are  forced  more  strongly  by 
the  winds,  the  waves  are  generally  steeper  and  hence  break  at  larger  scales  closer  to 
the  peak  in  the  wave  spectrum.  Because  breaker  momentum  increases  as  c1  the  total 
number  of  breakers  required  to  balance  the  wind  stress  decreases  with  decreasing 
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Figure  3.  Variation  of  the  peak  phase  speed  for  breaker  momentum  flux  cr  (solid  line)  and 
energy  flux  ce  (dashed  line)  with  wind  speed  and  wave  age.  The  phase  speeds  are  normalized 
by  the  maximum  phase  speed  in  the  wave-height  spectrum  cp. 


wave  age.  Meanwhile  as  winds  and  waves  approach  equilibrium  the  overall  steepness 
of  the  wave  field  relaxes  and  then  breaking  shifts  towards  the  smaller  (steeper)  waves. 
Again  because  of  the  c7  dependence  the  total  number  of  breakers  required  to  support 
the  wind  stress  in  the  equilibrium  range  must  then  increase,  i.e.  many  small  breakers 
dump  their  momentum  into  currents.  Larger  breakers  are  more  intermittent  but  their 
associated  vorticity  fields  are  more  vigorous  and  capable  of  interacting  with  the  vortex 
force.  The  interaction  between  breaking  and  vortex  force  thus  exhibits  a  wave  age 
dependence. 


4.  LES  model  of  the  OBL  with  wave  effects 

The  LES  model  we  use  to  examine  OBL  turbulence  is  a  conventional  one  (Moeng 
1984;  Sullivan,  McWilliams  &  Moeng  1996;  McWilliams  et  al.  1997)  but  is  extensively 
modified  to  account  for  surface  wave  effects:  it  is  based  on  the  incompressible 
Boussinesq  equations  with  a  single-point  second-moment  turbulent  kinetic  energy 
(TKE)  closure  subgrid-scale  parameterization  and  a  flat  top  (ocean)  surface.  The 
added  wave  effects  are  the  vortex  force  and  Lagrangian  mean  advection  associated 
with  Stokes  drift,  us',  and  a  wave-averaged  increment  to  the  pressure  that  arises 
through  conservative  wave-current  interaction  (McWilliams  et  al.  1997),  as  well  as 
additional  acceleration  and  energy  generation  due  to  non-conservative  wave  breaking 
(Sullivan  et  al.  2004).  The  governing  equations  for  momentum,  subgrid-scale  TKE, 
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and  density  (i.e.  the  same  as  for  any  material  scalar  concentration)  are  the  following: 
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(4.16) 


(4.1c) 


In  (4.1)  spatially  hltered  (or  resolved)  variables  are  denoted  with  an  overbar.  Here 
(u.w  =  Vxu)  are  the  resolved-scale  velocity  and  vorticity;  e  is  the  subgrid-scale  TKE; 
p  is  the  density;  and  the  generalized  pressure  held  is 


2  1  rr- 

3e+2 


+  nff)  ( Ui  +  uft') 


—  UiUi 


(4.2) 


Note  with  this  hux  form  of  the  momentum  equations,  n  reduces  to  ~p/p0  +  2c/ 3  in 
the  absence  of  a  wave  held.  The  effects  of  m  =  1  ,  ...,M  discrete  wave-breaking 
events  are  represented  by  a  resolved-scale  acceleration.  A",  and  a  subgrid-scale  TKE 
generation  rate,  Wm.  The  dots  in  the  TKE  equation  denote  all  the  other  terms 
(namely  advection,  production,  buoyancy,  diffusion  and  dissipation)  that  appear  in 
the  conventional  closure  formulation.  Other  variables  appearing  in  this  equation 
set  are:  the  Coriolis  frequency,  /  =  (0,0,/);  reference  density,  p0;  gravitational 
acceleration,  g;  and  subgrid-scale  momentum  and  density  huxes,  (t!;-,  r,p),  respectively. 
These  subgrid-scale  huxes  are  modeled  using  the  eddy  viscosity  prescription  described 
by  Moeng  (e.g.  Moeng  1984)  and  Sullivan  et  al.  (1994)  which  implies  that  the  principal 
feedback  of  e  on  u  and  p  occurs  through  subgrid-scale  mixing  with  eddy  viscosity 
and  diffusivities  (v,,  ly)  oc  e1/2. 

This  LES  model  (4.1)  includes  some  wave  effects  without  actually  resolving  the 
oscillatory  and  breaking  wave  motions  themselves.  By  otherwise  neglecting  wave 
dynamics,  there  is  an  implicit  assumption  that  the  wave  quantities  are  unaffected  by 
the  currents;  this  is  explicit  in  the  particular  rules  specihed  for  A,  W  and  uSl .  Further 
simplihcations  are  the  neglect  of  the  positive  buoyancy  caused  by  air  entrained  into 
bubbles  in  breaking  waves  (Lamarre  &  Melville  1991)  and  the  assumption  of  a 
linearized  equation  of  state  where  density  is  proportional  to  temperature  6.  Since 
the  wave  influences  on  deforming  the  oceanic  free  surface  are  averaged  out  in  this 
formulation,  the  wave  effects  on  pressure  do  not  enter  explicitly  into  the  model  since 
7i  is  calculated  by  solving  the  pressure-Poisson  equation  derived  from  taking  the 
divergence  of  the  momentum  equation  and  applying  the  incompressibility  condition. 
This  is  the  standard  method  of  determining  the  pressure  with  numerical  techniques 
based  on  fractional  step  methods  (Armfield  &  Street  1999;  Ferziger  &  Peric  2002). 


4.1.  Vortex  force 

Wave-averaged  influences  enter  as  Stokes  advection  in  the  TKE  and  scalar  equations 
and  as  a  generalized  vortex  force,  uSl  x  (fz+a i),  in  the  momentum  equation.  The  vortex 
force  has  a  significant  influence  on  the  turbulent  eddies  and  their  vertical  fluxes  in 
the  OBL  through  the  generation  of  Langmuir  circulations  (Craik  &  Leibovich  1976; 
McWilliams  et  al.  1997).  Wave  effects  also  directly  impact  subgrid-scale  energy  since 
its  balance  equation  includes  a  Stokes  production  term,  i.e.  subgrid-scale  stresses 
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working  on  the  vertical  gradient  of  the  Stokes  drift.  An  analogous  term  appears  in 
the  TKE  balance  for  resolved  scales  (McWilliams  et  al.  1997,  equation  5.1),  and  Stokes 
production  has  recently  been  included  in  some  one-dimensional,  Reynolds-averaged 
models  of  the  OBL,  (Kantha  &  Clayson  2004).  In  contrast  to  the  Reynolds-averaged 
models,  the  local  subgrid-scale  energy  transfer  from  Stokes  production  in  an  LES  can 
be  either  positive  or  negative  since  3 uSt /dz  has  a  given  orientation  and  sign  and  xn 
fluctuates  in  time  and  space  over  the  LES  grid.  Because  of  the  eddy-viscosity  closure 
for  vertical  momentum  flux,  the  long-time  impact  of  Stokes  production  for  e  depends 
on  the  current  and  Stokes  drift  shears;  i.e.  for  a  steady  uniform  uSt, 


(  du  3  w\ 
V  dz  +  dx  ) 


3  us'  / 

“37  +  V 


dv  dw 
3  z  +  3  y 


dvSt 

dz 


(4.3) 


where  (  )  indicates  an  ensemble  average.  The  magnitude  of  the  Stokes  drift  often 
exceeds  the  horizontal  current  speed  near  the  surface,  leading  to  substantial  TKE 
advection.  Finally,  the  coupling  of  uSr  and  e  is  implicit  because  of  the  influence  of 
waves  on  the  resolved  currents.  Stokes  advection  and  production  of  subgrid-scale  e, 
which  are  absent  in  our  previous  LES  (McWilliams  et  al.  1997),  are  expected  to  be 
more  important  in  the  current  simulations  that  account  for  a  full  wave  spectrum. 


4.2.  Resolved  and  subgrid-scale  breakers 

In  order  to  implement  our  LES  with  breakers  we  need  to  account  for  the  spatial 
filtering  associated  with  numerical  discretization.  Formally,  a  filtered  field  fix)  — *  fix) 
is  defined  by  the  operation, 


fix)  =  j  fix)Gix-x)dx, 


(4.4) 


where  G(x)  is  a  low-pass  spatial  filter  with  characteristic  length  scale  A ;  resolved  and 
subgrid-scale  motions  have  scales  /  >  A  and  /  <  A,  respectively. 

The  computational  grid  spacing  limits  our  ability  to  temporally  and  spatially 
resolve  the  breaker  wave  field.  Thus  we  introduce  a  minimum  breaker  phase  speed 
cA  =  igA/2n)1/2,  based  on  the  LES  filter  scale  A,  to  separate  resolved  and  unresolved 
breakers.  The  intermittent  space-time  properties  of  breakers  c  <  cA  cannot  be 
captured  on  the  LES  grid,  but  their  bulk  influence  must  still  be  accounted  for 
as  they  contribute  to  the  overall  fluxes  of  momentum  and  energy.  A  consistent  scale 
decomposition  of  the  momentum  flux  due  to  breaking  is  then 

r a  =  N  I  A  Jfic)P(c)dc  +  N  r  Jt{c)P{c)dc,  (4.5) 

J  0  J ca 

where  the  first  and  second  terms  on  the  right-hand-side  are  contributions  from 
subgrid-scale  and  resolved  breakers,  respectively.  Equation  (4.5)  follows  from 
combining  (3.4),  (3.6)  and  (3.7)  and  the  discussion  in  §3.3.  In  computational  practice, 
any  draws  made  from  the  PDF  with  phase  speed  c  <  cA  are  considered  subgrid-scale 
and  parameterized  as  a  constant  (viscous)  stress  x„  applied  at  the  water  surface.  Thus 
the  LES  representation  of  (4.5)  is 

rep 

xa  —  x0  N  I  Jt[c)  P(c)dc .  (4.6) 

J  CA 

The  subgrid-scale  contribution  to  momentum  flux  xa  is  known  when  the  grid  mesh 
is  specified,  i.e.  x0  is  calculated  prior  to  integration  of  the  governing  equations.  Note 
(4.5)  behaves  properly  with  changing  grid  spacing.  In  the  limit  where  A  is  large  all 
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breaking  becomes  subgrid-scale  and  then  the  surface  forcing  reduces  to  conventional 
constant  stress. 

By  analogy  with  momentum  flux  the  decomposition  of  the  wave  (breaker)  energy 
flux  into  small-  and  large-scale  pieces  is 

rCp 

Ea  =  E0  +  N  c-Jt(c)P(c)dc,  (4.7) 

J  Ca 

where 

£.„  =  N  [  *  c-J?(c)P(c)dc.  (4.8) 

Jo 

The  energy  flux  £„  from  small-scale  breakers  c  <  cA  is  tightly  confined  to  the  water 
surface,  while  the  energy  flux  from  larger-scale  breakers  c  >  cA  is  naturally  distributed 
over  a  finite  horizontal  extent  and  vertical  depth  spanning  multiple  gridpoints. 

In  order  to  expose  the  breaker  energetics  in  our  LES  model  we  construct 
kinetic  energy  equations  for  the  total,  resolved  and  subgrid-scale  components.  The 
contribution  of  breakers  to  the  evolution  of  filtered  total  kinetic  energy  per  unit  of 
water  mass  E  =  WTuj/2, 


dE 

~dt 


dui 

~dt~ 


+  W;  At  , 


(4.9) 


depends  on  the  (unknown)  correlation  between  total  velocity  and  breaker  impulses. 
Elere  dots  denote  all  other  terms  in  the  kinetic  energy  equation,  i.e.  advection, 
production,  buoyancy,  diffusion  and  dissipation  (Moeng  1984).  Combining  (4.9)  with 
the  transport  equation  for  resolved  kinetic  energy  Er  =  uj TTJ/2, 


dEr 

~dT 


Ui 


duj 

Jdi 


+  Ui  A[ , 


(4.10) 


leads  to  the  TKE  equation  for  subgrid-scale  energy  e  =  (m,  u,  —  u ,•  «,■)/ 2 


de  _  dE  dEr 
dt  d t  dt 

In  (4.11),  the  subgrid-scale  work  done  by  breaking, 

W  =  Ui  Aj  —  ujAi , 


(4.11) 

(4.12) 


is  unknown  and  must  be  modeled  consistently  in  terms  of  the  resolved  field  to  match 
the  input  energy  flux  £a.  Equating  the  volume  work  done  by  breakers  in  the  water 
column  to  the  input  energy  flux  from  winds  to  waves  requires 

—  =  4  f  u‘A,dr  =  1  f  (ujAi  +  W)dr.  (4.13) 

Po  '  ,/  i‘  '  J  i‘ 

Inspection  of  (4.13),  (4.7),  and  (4.8)  guides  our  model  design  for  W.  We  adopt  the 
two-part  model, 

W(x,  t)  =  —  ^H(z)  +  V  kw  dn  -  A" ,  (4.14) 

po  dz  ^ 

m 

where  the  Heaviside  function  is  defined  as  H(z  ^  0)  =  1,  H(z  <  0)  =  0.  The  first  term 
in  (4.14)  is  the  gradient  of  the  small-scale  energy  flux  from  the  wave  field  at  the  water 
surface  (z,  =  0)  and  can  be  explicitly  computed  from  the  discrete  event  breaker  model 
and  the  PDF  of  breaking;  it  varies  with  wind  speed,  wave  age  and  LES  filter  scale 
A.  The  second  term  in  (4.14)  is  the  work  done  by  the  sum  of  m  individual  large-scale 
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breakers  with  c  >  cA ;  it  distributes  breaker  work  over  a  finite  depth  and  mimics  the 
intermittent  space-time  character  of  the  momentum  impulses.  In  principle,  given  S'a 
and  £0  the  work  constant  kw  is  chosen  to  satisfy  (4.13).  However,  because  of  the 
dependence  on  the  resolved  field  in  (4.13),  kw  can  only  be  found  iteratively,  i.e.  by 
first  choosing  a  trial  kw,  running  the  LES  code,  and  post-checking  the  results  against 
S’a.  Then  kw  is  adjusted  accordingly  to  satisfy  (4.13);  we  find  kw  ~  1.3. 

It  is  interesting  to  compare  our  model  for  breaker  work  (4.14)  with  the  breaker 
energy  flux  model  first  proposed  by  Craig  &  Banner  (1994)  (the  latter  is  used  in  the 
LES  described  by  Noh  et  al.  (2004)).  Craig  &  Banner  (1994)  collapse  the  influence  of 
breaking  onto  the  water  surface,  with  no  depth  dependence,  and  replace  the  temporal 
and  spatial  variability  of  breaking  by  an  ensemble  average  inherent  in  second-order 
closure  modeling.  The  first  term  of  (4.14)  is  similar  to  the  Craig  &  Banner  (1994) 
model  but  only  has  contributions  from  the  subgrid-scale  piece  of  the  wave  spectrum, 
i.e.  breakers  with  c  <  cA.  This  is  consistent  with  the  LES  decomposition  into  resolved 
and  subgrid-scale  fields.  For  a  fixed  wave  state,  as  A  — >  oo  the  work  done  by  breaking 
tends  to  the  ensemble  average.  Unlike  Craig  &  Banner  (1994)  £0  varies  with  wave 
age  because  of  the  dependence  of  the  PDF  on  wave  state. 


5.  Simulations  with  wave  effects 

In  order  to  evaluate  the  impacts  and  interactions  between  wave  breaking,  vortex 
force  and  turbulence  we  incorporated  the  modeling  ideas  and  equations  of  §  2-4  into 
a  large-eddy  simulation  code  for  the  OBL.  Similarly  to  McWilliams  et  al.  (1997)  the 
boundary  conditions  are  periodic  in  the  lateral  (x—y)  directions,  no  flow  at  the  lower 
boundary  with  a  radiation  condition  allowing  gravity  waves  to  escape  (Klemp  & 
Durran  1983;  McWilliams  et  al.  1997).  The  imposed  surface  stress  depends  on  wind 
speed  and  the  PDF  of  breaking  as  discussed  previously.  The  LES  algorithm  with 
stochastic  breakers  is  sufficiently  novel  (and  complex)  that  a  brief  description  of  the 
code  is  provided  in  the  Appendix.  Further  algorithmic  details,  numerical  method  and 
the  subgrid-scale  model  used  to  solve  the  equation  set  (without  wave  breaking)  are 
fully  presented  in  McWilliams  et  al.  (1997),  Sullivan  et  al.  (1996),  Sullivan,  McWilliams 
&  Moeng  (1994)  and  Moeng  (1984).  A  direct  numerical  simulation  code  for  the  OBL 
with  single-scale  breakers  is  outlined  in  Sullivan  et  al.  (2004). 

An  exhaustive  exploration  of  the  large  parameter  space  encompassing  variations 
in  wind  speed,  wave  age,  buoyancy,  varying  combinations  of  wave  effects  and  other 
physical  processes  is  not  possible  because  of  finite  computational  resources.  The  scope 
of  the  current  study  focuses  on  OBL  dynamics  when  wave  influences  are  important 
and  in  particular  when  intermittent  breaking  is  a  dominant  surface-layer  process.  For 
the  present  investigation  we  choose  a  moderately  high  wind  speed  Ua  =  15  m  s-1. 
The  suite  of  experiments  compares  four  types  of  simulations:  a  baseline  simulation 
driven  by  uniform  (constant)  surface  stress  with  no  wave  effects;  the  conventional 
posing  plus  wave-averaged  Stokes  drift  terms  (McWilliams  et  al.  1997);  a  posing 
with  stochastic  forcing  representing  breaking  waves;  and,  a  posing  with  stochastic 
forcing  plus  Stokes  drift  terms.  In  addition,  solution  sensitivity  tests  are  performed 
to  examine  the  consequences  of  modestly  varying  wave  age  (as  contained  in  the 
Terray  parameter),  an  increase  in  wind  speed  to  Ua  =  30  m  s_1  and  a  decrease  in 
the  Stokes  drift  profile  for  a  fetch-limited  wave  spectrum.  These  explorations  are 
intended  to  illustrate  the  solution  robustness  and  are  mild  excursions  into  the  regime 
of  disequilibrium  winds  and  waves.  A  summary  of  the  simulations  and  the  naming 
conventions  used  to  identify  the  various  cases  is  given  in  table  1.  A  wave  age  range  is 
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Run  name 

Forcing 

Stokes  drift 

Cp/u*a 

gt 

Ua  (m  s  1 

)  u,„  (ms1 

)  h  (m) 

U 

uniform  stress 
no  wave  effects 

NA 

NA 

NA 

15 

0.0187 

-48.6 

U  +  St 

uniform  stress 

Yes 

NA 

NA 

15 

0.0187 

-51.3 

B 1 

breaking 

No 

>30 

3.9 

15 

0.0187 

-46.6 

B 1  +  St 

breaking 

Yes 

>30 

3.9 

15 

0.0185 

-48.6 

B  2 

breaking 

No 

23  +  3 

5.5 

15 

0.0184 

-49.4 

B2+St 

breaking 

Yes 

23  +  3 

5.5 

15 

0.0187 

-46.3 

B  3 

breaking 

No 

19  +  3 

7.4 

15 

0.0187 

-52.1 

S3  +  St 

breaking 

Yes 

19  +  3 

7.4 

15 

0.0195 

-48.4 

•[S3  +  FSt 

breaking 

Yes 

19  +  3 

7.4 

15 

0.0194 

-46.7 

HI  +  St 

breaking 

Yes 

>30 

3.9 

30 

0.0451 

-71.7 

Table  1.  Simulation  properties 


f  The  Stokes  drift  in  this  simulation  is  based  on  a  fetch-limited  wave-height  spectrum  as  discussed 

in  §7. 


supplied  in  table  1  to  indicate  the  uncertainty  in  the  estimate  of  g,  based  on  the  data 
of  Terray  et  al.  (1996,  figure  8). 

The  setup  of  the  simulations  is  similar  to  McWilliams  et  al.  (1997)  except  we  use 
larger  domains  and  finer  grid  resolutions.  The  initial  state  of  the  OBL  is  neutral 
stratification  from  the  surface  to  a  mixed  layer  depth  h  =  —32  m.  Below  z  <  h  the 
therm ocline  is  stably  stratified  at  a  rate  of  0.05  K  m_1.  In  order  to  examine  scalar 
mixing  a  small  heat  flux  Q,  =  5  x  10~7  K  m  s_1  is  imposed  at  the  surface.  In  all 
simulations  the  Coriolis  parameter  /  =  10~4  s~3.  For  the  wind  speed  Ua  =  15  m  s_1 
and  assuming  a  density  ratio  p0/pa  =  103  the  atmosphere  and  ocean  friction  velocities 
are  u,a  =  0.591  m  s-1  and  u,0  =  0.0187  m  s_1,  respectively.  The  computational 
domain  ( Lx ,  Ly,  Lz)  =  (300,  300,  —110)  m  is  discretized  with  (300,  300,  128)  gridpoints 
(horizontal  grid  spacing  Ax  =  Ay  =  1.0  m).  At  the  higher  wind  speed  Ua  =  30  m  s_1, 
a  larger  domain  is  used,  (750,  750,  —170)  m,  and  is  discretized  with  (500,  500,  128) 
gridpoints  (Ax  =  Ay  =  1.5  m).  In  this  case  the  rate  of  stable  stratification  below 
the  initial  thermocline  depth  is  doubled  0.10  K  m_1.  Vertical  meshes  are  generated 
using  constant  algebraic  stretching,  i.e.  with  the  ratio  of  any  two  vertical  cells  K  = 
Azk+i/Azk  held  constant.  Stretching  factors  K  <  1.012  yield  flexible  meshes  capable 
of  concentrating  the  grid  near  the  surface,  but  extending  deep  into  the  water.  This 
small  value  of  K  ensures  that  the  mesh  varies  smoothly  in  the  vertical  direction. 
Thus  the  ratio  of  horizontal  and  vertical  spacing  at  the  first  z-level  below  the  water 
surface  Ax/Azi  <  2.6.  LES  with  this  horizontal-vertical  mesh  ratio  yields  good 
simulations  of  atmospheric  convection  (Moeng  &  Wyngaard  1989)  and  this  mesh 
ratio  is  well  within  the  acceptable  limits  suggested  by  Scotti,  Meneveau  &  Lilly 
(1993).  The  influence  of  the  subgrid-scale  (SGS)  motions  is  not  significant  as  the  SGS 
fluxes  are  small  over  the  entire  computational  domain,  especially  in  simulations  with 
wave  breaking.  To  arrive  at  a  (statistical)  steady  state  all  simulations  are  run  for 
more  than  30  large  eddy  turnover  times  Te  =  —h/u,0  or  more  than  25  physical  hours 
requiring  [90000—240000]  computational  steps.  As  is  customary  practice,  statistics 
are  generated  by  combining  spatial  x  — y  and  temporal  averaging;  these  averages  are 
indicated  by  (  ).  The  mixed-layer  depth  is  a  function  of  time  and  also  varies  with 
spatial  position.  At  any  t,  we  deduce  an  estimate  of  the  mean  mixed-layer  depth  h 
using  the  maximum  scalar  gradient  method  described  by  Sullivan  et  cd.  (1998).  This 
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method  is  able  to  track  the  average  movement  of  the  turbulent  layer  and  provides  a 
good  estimate  of  the  average  entrainment  velocity  across  the  thermocline. 

With  wave  effects  the  vertical  profile  of  the  Stokes  drift  is  numerically  evaluated 
for  a  full  wave  spectrum  as  described  in  §2.1.  Based  on  the  Stokes  pro  hie  at  a  depth 
z  =  1  m,  the  turbulent  Langmuir  number  La,  =  \J u»0/uSl  «  0.3,  and  thus  wave- 
current  interactions  are  important  (McWilliams  et  al.  1997)  in  the  chosen  wind-speed 
regime.  The  partitioning  of  the  atmospheric  momentum  flux  (2.5)  into  resolved  and 
subgrid  breaker  pieces  depends  on  the  grid  resolution,  wind  speed  and  the  breaker 
PDF  as  discussed  in  §§3.3  and  4.2.  For  the  grid  resolution  used,  winds  Ua  =  15  m  s_1, 
and  near  equilibrium  winds  and  waves  g,  «  3.9,  the  ratio  of  the  subgrid-scale  uniform 
stress  to  the  total  wind  stress  |r0|/|ra|  «  0.1.  For  higher  wind  speeds  and  larger  values 
of  g,  this  ratio  rapidly  decreases  to  values  less  than  1%.  Hence  for  the  majority  of 
the  simulations  intermittent  breaking  is  responsible  for  transmitting  90-100%  of  the 
atmospheric  stress  to  the  water  column  and  is  the  main  mechanism  for  mean  current 
generation  in  the  simulations.  Below  a  depth  z  <  —4  m  the  SGS  contribution  to  the 
momentum  fluxes  is  less  than  2%  of  the  surface  stress  n;0. 


6.  Results 

The  observed  response  of  the  OBL  to  the  vortex  force  and  wave  breaking  suggests 
a  connection  between  the  surface  wave  conditions  and  vertical  mixing  in  the  OBL.  We 
find  that  the  mean  currents,  turbulence  variances  (and  TKE),  turbulence  dissipation, 
scalar  and  momentum  fluxes  and  entrainment  at  the  thermocline  all  exhibit  varying 
degrees  of  sensitivity  to  the  surface  wave  field.  In  general,  simulations  with  explicit 
wave  influences  differ  from  their  counterparts  driven  by  constant  stress  and  no  wave 
effects.  An  important  parameter  is  the  PDF  of  breaking  and  its  dependence  on  wave 
state.  For  a  given  wind  speed  P(c)  shifts  towards  larger-scale  breaking  as  wave  age 
cp/u,a  decreases  from  30  to  19.  At  the  same  time  the  number  of  breakers  in  the 
computational  domain  decreases,  consistent  with  the  need  to  obey  momentum  and 
energy  conservation.  The  shift  towards  larger-scale  breaking  promotes  interactions 
between  breaking  turbulence  and  vortex  force. 

6.1.  Mean  current  and  momentum-flux  profiles 
In  order  to  examine  the  momentum  balance  in  the  presence  of  waves  we  form  the 
ensemble-average  integral  momentum  budget  from  the  LES  equations  (4.1a).  For  the 
conditions  of  the  simulations,  i.e.  horizontally  homogeneous  flow,  no  flow  as  z  — *  — oo, 
aligned  surface  winds  and  waves  with  uSt  =  [n5'(z),  0,  0],  the  momentum  budget  for 
statistically  steady  flow  is 

(urn)  +  (ri3)  -  f  (A)  dz  =  f  f  (v)  d Z  (6.1a) 

J  —CO  J  —oo 

(vw)  +  (r23>  =  —  /  f  (u  +us,)d z,  (6.1  b) 

J  —CO 

where  (  )  denotes  an  ensemble  average.  This  expression  illustrates  how  resolved 
and  subgrid-scale  turbulence  fluxes,  waves  and  mean  currents  contribute  to  the  bulk 
momentum  balance.  It  contains  explicit  wave  influences  in  the  form  of  Ekman-Stokes 
transport  f  us,dz  (McWilliams  et  al.  1997)  and  a  new  contribution  from  vertically 
distributed  breakers.  Note  in  the  y-budget  equation  the  locally  varying  vortex  force 
us'a>i,  which  plays  an  important  role  in  generating  convergence  lines,  does  not  appear 
since  ( us'co 3)  =  0.  The  vertical  profile  of  average  breaker  momentum  flux  —f(A)dz 
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Figure  4.  The  depth  variation  of  the  ensemble-average  breaker  momentum  flux  for  varying 
wave  age  and  wind  speed:  (Ua  ,cp/u,a )  =(15  m  s  30),  •;  (15  m  s  23),  o;  (15  m  s  *,  19), 
□  ;  (30  m  s_1,  30),  A.  The  breaker  momentum  flux  is  normalized  by  the  water  friction  velocity 


obtained  from  the  simulations  is  shown  in  figure  4.  As  expected  its  surface  value  varies 
with  wave  age  and  wind  speed,  as  discussed  in  §  3.  Intermittent  breakers  account  for 
more  than  90%  of  the  total  momentum  transport  to  the  underlying  currents  for  the 
forcing  conditions  and  grid  resolution  of  the  simulations.  The  breaker  momentum 
flux  decays  rapidly  with  z  because  of  the  exponential  decay  of  the  vertical  shape 
function  (Sullivan  et  al.  2004,  (3.3d))  and  the  small  value  of  the  depth  penetration 
constant  /  =  0.2;  for  example  the  breaker  momentum  flux  with  Ua  =  15 ms-1  and 
gt  =  3.9  becomes  negligible  below  z  *  —  lm  and  at  higher  winds,  where  the  PDF 
includes  larger  breakers,  extends  down  to  about  z  «  —5  m.  Notice  for  intermediate 
wave  age  the  breaker  contribution  at  15 ms-1  moves  towards  the  variation  at  30 ms-1 
consistent  with  the  PDF  variations  shown  in  figure  1. 

In  order  to  compare  fairly  average  currents  from  different  simulations,  which  have 
slightly  different  integration  periods,  the  mean  currents  are  corrected  for  the  presence 
of  inertial  oscillations.  Given  the  Coriolis  parameter,  the  amplitude  and  phase  of  the 
inertial  oscillation  is  determined  by  a  least-squares  curve-fitting  procedure  applied 
over  the  entire  time  period  of  the  simulations  (see  Lin  et  al.  1996).  The  deduced  low- 
frequency  inertial  oscillation  is  then  subtracted  from  each  current  component.  Vertical 
profiles  of  the  mean  currents,  shown  in  figure  5,  illustrate  the  impact  of  the  various 
wave  processes  under  conditions  approaching  a  fully-developed  sea.  The  most  striking 
feature  of  the  results  is  a  clear  grouping  based  on  the  presence  or  absence  of  Stokes 
drift.  For  wind-wave  equilibrium  conditions,  breaking  does  not  significantly  modify 
the  mean  direction  and  magnitude  of  the  average  currents.  This  is  an  interesting  result 
given  the  radical  difference  in  the  level  of  intermittent  forcing  between  simulations 
with  and  without  breaking.  We  emphasize  that  in  cases  (51,  51  +  St)  mean  currents 
are  generated  solely  by  spatially  distributed  random  impulses  near  the  water  surface, 
while  in  cases  (U,  U  +  St)  a  conventional  uniform  surface  stress  boundary  condition 
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( v)/u » 


Figure  5.  The  sum  of  the  mean  current  and  Stokes  drift  profiles  for  simulations  with:  no 
wave  effects,  dash-dot  line;  Stokes  drift  only,  ►;  breaking  only,  solid  line;  Stokes  drift  plus 
breaking,  •.  The  inset  shows  just  the  Eulerian  mean  current  profile  (u)/u,0  for  simulations 
with  Stokes  drift  only  (►)  and  Stokes  drift  plus  breaking  (•).  The  winds  and  waves  are  near 
full  equilibrium  with  ( cp/iua ,  g, )  =  (30,  3.9). 


is  applied  at  z  =  0.  Under  the  assumption  that  breaking  is  the  main  path  to  current 
generation,  these  results  imply  the  cumulative  effect  of  intermittent  breaking  is  to 
generate  mean  currents  typical  of  an  OBL  forced  by  constant  surface  stress.  Closer 
inspection  of  the  results  shows  the  mean  surface  current  obtained  with  breaking  only 
(no  vortex  force)  weakens  as  wave  age  decreases.  As  cp/u,a  diminishes  the  probability 
of  large-scale  breaking  increases,  which  induces  a  substantial  increase  in  near-surface 
eddy  viscosity  v,  accompanied  by  a  reduction  in  the  current  shear.  A  similar  effect 
was  also  observed  in  our  idealized  DNS  (Sullivan  et  al.  2004)  but  is  weaker  in  the 
present  LES  which  includes  a  spectrum  of  breaking  waves. 

The  presence  of  Stokes  drift,  however,  has  a  dramatic  impact  on  the  mean  currents 
and  plays  a  dominant  role  even  as  the  level  of  breaking  increases.  The  present  results 
are  at  least  qualitatively  consistent  with  the  previous  LES  of  McWilliams  et  al.  (1997), 
who  applied  a  Stokes  profile  based  on  a  monochromatic  wave  field  and  drove  the 
OBL  with  uniform  surface  stress.  Potent  Langmuir  cells,  generated  by  the  wave 
field,  promote  efficient  vertical  transport  leading  to  well-mixed  ( u ,  u)-current  profiles 
throughout  the  bulk  of  the  mixed  layer  with  a  noticeable  turn  to  the  right  as  z  — *  h. 
Langmuir  cells  also  indirectly  impact  entrainment  of  cool  water  as  they  enhance  the 
current  shear  in  the  region  z  «  h. 

The  coherent  structures  generated  by  the  vortex  force  increase  the  efficiency  of  the 
momentum  transport  as  observed  in  figure  6.  For  a  given  z,  the  total  turbulent  fluxes 
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Figure  6.  Vertical  profiles  of  the  mean  total  (resolved  plus  SGS)  vertical  momentum  flux 
(i u!w )  for  simulations  with:  no  wave  effects,  dash-dot  line;  Stokes  drift  only,  ►;  breaking  only, 
solid  line;  Stokes  drift  plus  breaking,  •. 

( u'w )  and  ( v'w )  are  enhanced  compared  to  the  situation  with  no  wave  influences 
and  this  enhancement  is  only  slightly  modified  by  the  presence  of  significant  wave 
breaking  for  fully  developed  waves.  The  profile  (figure  7)  of  mean  eddy  viscosity  for 
momentum  Km,  deduced  from 

(u'w)  =  -Kmd(u)/dz,  (6.2) 

illustrates  the  increased  mixing  efficiency  generated  by  Langmuir  cells.  The 
enhancement  compared  to  cases  with  no  vortex  force  is  primarily  a  reflection  of 
the  decrease  in  mean  current  shear  shown  in  figure  5.  The  sign  change  of  Km  in 
the  interval  —0.2  <  z/\h\  <  0  is  a  consequence  of  subtle  slope  changes  in  the  (v) 
and  (if)  currents  near  z/\h\  «  (—0.2, —0.05),  respectively.  Vertical  oscillations  in  the 
K,„  profile  in  the  case  with  vortex  force  appear  to  be  smoothed  by  the  presence  of 
stochastic  breaking.  The  negative  values  of  the  mixing  coefficient  are  evidence  of 
non-local  vertical  transport  and  clearly  show  the  inadequacy  of  a  mean  eddy  viscosity 
assumption  in  the  presence  of  Langmuir  cells. 

6.2.  Variance,  TKE  and  dissipation  profiles 

Inspection  of  the  TKE  and  turbulence  variance  profiles  (figures  8  and  9)  highlights 
the  competition  and  tradeoffs  between  uniform  and  stochastic  forcing,  and  Stokes 
drift  at  wind  speed  Ua  =  15 ms-1.  Similar  to  the  mean  current  profiles,  the  results 
for  the  resolved  variances  can  be  grouped  depending  on  the  presence  or  absence  of 
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Figure  7.  Profile  of  eddy  viscosity  for  momentum  Km  =  — ( u'w )  •  d(«)/dz  /  |d(S)/dz|2  for 
simulations  with:  no  wave  effects,  dash-dot  line;  Stokes  drift  only,  ►;  breaking  only,  solid  line; 
Stokes  drift  plus  breaking,  •. 

vortex  force.  First,  comparing  flows  driven  by  breaking  or  uniform  stress  (no  Stokes 
drift),  we  notice  that  the  influence  of  breaking  is  confined  to  the  near-surface  waters 
—0.1  <  z/\h\  <  0,  a  depth  of  5  m  or  less  that  is  expected  based  on  the  vertical 
distribution  of  breaker  impulses  (figure  4).  In  the  LES,  stochastic  breaking  lowers  the 
u -variance  similarly  to  the  trend  observed  in  our  DNS,  which  had  relatively  large 
breakers  of  a  single  scale.  The  biggest  impact  of  the  wave  field  on  the  turbulence 
results  from  the  vortex  force.  The  near-surface  reduction  in  the  u  -variance,  strong 
enhancement  of  the  u-variance  and  the  significant  increase  in  the  w-variance  below 
the  water  surface  observed  here  are  also  present  in  the  LES  of  McWilliams  et  al. 
(1997)  obtained  with  a  monochromatic  wavefield,.  These  changes  in  the  turbulence 
variances  are  attributable  to  the  presence  of  organized  Langmuir  cells  which  persist 
in  the  presence  of  breaking.  Closer  inspection  of  runs  B 1  +  St  and  U  +  St  indicates  a 
reduction  in  the  v -variance  in  the  presence  of  breaking  near  the  water  surface.  This 
trend  becomes  more  apparent  as  the  dominant  scale  of  breaking  increases,  i.e.  for 
less  developed  waves  cp/u,a  <  30  (results  not  shown).  The  reduction  in  the  lateral 
variance  implies  less  coherent  Langmuir  cells;  this  finding  is  supported  by  the  flow 
visualization  presented  later.  The  modulation  of  the  strength  and  organization  of 
Langmuir  cells  by  large-scale  breaking  waves  in  the  LES  is  an  interesting  result  as 
previous  LES  show  strong  intensification  and  reduction  in  scale  of  the  Langmuir  cells 
as  the  LES  grid  resolution  is  refined  (McWilliams  et  al.  1997). 

Breaking  and  the  vortex  force  elevate  the  TKE  near  the  water  surface  compared  to 
simulations  without  wave  influences;  the  total  (resolved  plus  SGS)  energy  normalized 
by  ula  increases  from  6  to  more  than  20  with  a  large  percentage  of  the  increase  due  to 
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Figure  8.  TKE  profiles  close  to  the  water  surface  —0.2  <  z/\h\  <  0  for  simulations  with-no 
wave  effects,  dash-dot  line;  Stokes  drift  only,  ►;  breaking  only,  solid  line;  Stokes  drift 
plus  breaking,  •.  The  wave  age  is  cp/u,a  =  30.  Panel  (a)  total  (resolved  plus  SGS)  and 
\b)  subgrid-scale  energy. 


Figure  9.  Variance  profiles  of  the  resolved  velocity  components  for  simulations  with:  no  wave 
effects,  dash-dot  line;  Stokes  drift  only,  ►;  breaking  only,  solid  line;  Stokes  drift  plus  breaking 
for  winds  and  waves  near  equilibrium,  •. 
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Figure  10.  Profile  of  normalized  dissipation  e\h\/ula  for  simulations  with:  no  wave  effects, 
dash-dot  line;  Stokes  drift  only,  ►;  breaking  only,  solid  line;  Stokes  drift  plus  breaking,  •; 
fully-developed  wind  waves.  Dissipation  is  estimated  from  the  subgrid-scale  parameterization 
e  =  Cfe3^2/2\  with  Ce  =  0.93.  The  heavy  black  line  is  the  rough  wall  layer  scaling  e  =  w»o/0.4|z|. 
The  inset  is  a  blow-up  near  the  surface. 

the  enhancement  of  the  SGS  energy  by  breaking.  Note  in  the  region  0  >  z/\h\  >  —0.05 
breaking  is  overwhelmingly  the  dominant  source  of  SGS  energy.  Below  a  depth 
z/\h\  <  —0.1  the  simulations  again  segregate  into  two  families  based  on  the  presence 
or  absence  of  the  vortex  force  as  Langmuir  cells  are  depth  filling  and  enhance  the 
TKE  over  the  bulk  of  the  mixed  layer. 

Turbulent  dissipation  is  a  metric,  produced  by  observational  datasets,  often  used 
to  infer  the  impact  of  waves  on  ocean  currents,  e.g.  Terray  et  al.  (1996),  Drennan 
et  al.  (1996),  Terray,  Drennan  &  Donelan  (1999).  These  field  measurements  report 
elevated  dissipation  by  1-2  orders  of  magnitude  compared  to  (atmospheric)  wall 
scaling  e  ~  l/z  with  the  magnitude  dependent  on  wave  state.  This  increase  is  most 
often  attributed  to  breaking  waves.  Previous  LES  with  the  vortex  force  driven  by 
uniform  stress  also  find  enhanced  dissipation  near  the  water  surface,  and  thus  it 
is  interesting  to  examine  the  relative  contributions  of  Langmuir  cells  and  breaking 
waves  to  dissipation  in  the  present  solutions.  The  turbulent  dissipation  from  cases 
with  uniform  stress  and  different  combinations  of  vortex  force  and  wave  breaking 
are  compared  in  figures  10  and  11.  Here  we  use  the  traditional  SGS  parameterization 
e  =  Cee3/2 / A  (Moeng  1984)  as  a  measure  of  net  dissipation  in  our  simulations.  The 
LES  driven  by  uniform  stress  and  no  wave  effects  reasonably  reproduces  the  rough 
wall  estimate  e  =  u10/k\z\,  where  k  =  0.4. 

In  our  modeling,  breaker  work  is  strongly  intermittent  in  space  and  time  with 
its  magnitude  dependent  on  the  phase  speed  c  of  any  particular  breaker.  As  a 


Ocean  boundary  layers  with  vortex  force  and  stochastic  breaking 


431 


Figure  11.  Near-surface  normalized  dissipation  ep0h/Sa  for  simulations  with  breaking  and 
Stokes  drift  for  varying  wave  age;  cp/u,a  =  30,  dash-dot  line;  cp/u,a  =  23,  dashed  line;  and 
c p / u*a  =  19,  long-dashed  line.  The  measurements  of  Terray  et  al.  (1996)  and  Drennan  et  al. 
(1996)  are  indicated  by  •  and  A,  respectively.  In  these  data  the  length  scale  It,  =  cp/g,  where 
cp  is  the  peak  in  the  measured  wave  height  spectrum  for  wind  waves.  In  the  LES,  It,  =  c|/g, 
where  cE  is  the  phase  speed  of  the  peak  in  the  breaker  energy  flux  spectrum. 


consequence  of  (4.1  b),  all  SGS  variables  including  e  inherit  these  dependencies  and 
considerable  averaging  is  required  to  obtain  reliable  dissipation  estimates.  Below  a 
depth  z/\h\  <  —0.2  the  dissipation  profiles  for  all  cases  are  similar  in  shape,  and 
the  results  with  the  vortex  force  are  only  modestly  elevated  compared  to  cases 
with  no  Stokes  drift.  However,  near  the  water  surface  the  magnitude  and  vertical 
distribution  of  dissipation  depend  strongly  on  breaking  and  the  vortex  force.  With 
breaking  and  Stokes  drift  the  dissipation  close  to  the  water  surface  is  more  than 
(80,2)  times  larger  than  values  obtained  from  simulations  with  no  wave  effects  and 
Stokes  drift  only,  respectively.  Enhanced  dissipation  is  primarily  a  consequence  of 
breaking,  which  is  expected  based  on  the  high  values  of  TKE  shown  in  figure  8. 
Figure  11  compares  dissipation  estimates  from  LES  with  varying  wave  age  in  the 
non-dimensional  variables  (z/h)  and  (e4p0/£a)  proposed  by  Terray  et  al.  (1996)  and 
Drennan  et  al.  (1996),  where  4  is  a  characteristic  length  scale  of  the  wave  field.  In  the 
field  observations  e  is  obtained  from  velocity  spectra  in  the  inertial  subrange,  and  the 
majority  of  the  measurements  were  obtained  under  developing  seas  (wave  age  <  8) 
with  significant  wave  heights  Hs  <  0.5  m.  The  data  of  Drennan  et  cd.  (1996)  are  more 
developed  cp/u,a  <  23  with  Hs  ~  0.88—2  m.  In  terms  of  significant  wave  height,  all 
the  data  are  well  below  the  estimate  for  fully  developed  waves  discussed  in  §2.1  (Alves 
et  cd.  2003).  Drennan  et  al.  (1996)  note  the  choice  of  length  scale  in  the  dissipation 
scaling  law  is  not  unique;  they  plot  their  results  using  4  =  Hs  and  4  =  c2p/g,  where 
cp  is  the  peak  in  the  wave-height  spectrum  excluding  swell.  The  statistical  variations 
in  the  observations  are  large,  the  range  of  conditions  is  small  and  an  optimum  length 
scale  is  not  obvious  from  the  data.  These  uncertainties  make  a  direct  comparison  with 
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the  present  LES  results  obtained  at  higher  wind  speeds  difficult.  However,  to  illustrate 
the  trends  we  show  the  LES  results  normalized  using  a  length  scale  lb  =  c2E/g,  where 
cE  is  the  phase  speed  of  the  peak  in  the  breaker  energy-flux  spectrum  (see  figure  3). 
The  qualitative  comparison  between  the  observations  and  the  LES  with  breakers  and 
Stokes  drift  is  good.  As  the  wave  age  decreases  from  wind-wave  equilibrium  the  LES 
results  tend  to  limiting  values,  independent  of  cp/u,a,  in  agreement  with  the  Terray 
et  al.  (1996)  scaling  arguments.  A  wider  range  of  observational  datasets  is  needed  to 
test  the  LES  predictions. 


6.3.  Flow  structures 

The  statistical  moments  discussed  in  the  previous  sections  are  evidence  of  the  influence 
of  surface  waves  on  the  bulk  properties  of  the  mixed  layer.  Extensive  visualization  of 
the  various  flows  is  next  discussed  to  elucidate  further  the  spatial  and  temporal  effects 
of  vortex  force  and  breaking  on  mixed-layer  dynamics.  Here  the  vertical  velocity  field 
is  used  as  a  marker  to  identify  coherent  structures  in  the  flow.  Figure  12  compares 
w-contours  near  the  water  surface  for  flows  with  different  combinations  of  uniform 
stress  or  breaking  and  vortex  force.  The  baseline  simulation  with  no  breaking  and 
no  vortex  force  shows  that  randomly  distributed  turbulent  eddies  dominate  the  near¬ 
surface  motions,  as  is  typical  of  a  flat  plate  boundary  layer.  Intermittent  breaking  alone 
modifies  this  pattern  in  a  manner  consistent  with  our  previous  DNS:  each  breaker 
generates  a  forward  and  downward  impulse  and  a  weaker  positive  return  flow.  The 
breaker  flow  structures  scale  with  the  phase  speed  c  and  the  number  of  events  is 
consistent  with  the  PDF  of  breaking.  The  relatively  benign  patterns  in  figures  12(a) 
and  (c)  are  strongly  modulated  by  the  vortex  force.  Streamwise  elongated  patterns 
appear  in  figures  12(h)  and  ( d )  and  the  downwelling  and  upwelling  lines  reflect  the 
formation  of  Langmuir  cells.  Closer  inspection  of  the  patterns  reveals  streamwise 
mergers  at  forward-looking  Y-junctions,  a  consequence  of  the  merger  of  positive  and 
negative  signed  streamwise  vortices  from  neighboring  cells.  The  overall  strength  of  the 
downwelling  and  upwelling  lines  is  consistent  with  the  high  levels  of  vertical-velocity 
variance  shown  in  figure  9.  The  life  cycles  of  Langmuir  cells  for  mixed  layers  driven 
by  uniform  stress  are  fully  discussed  by  McWilliams  et  al.  (1997). 

There  is  a  subtle  hint  in  the  images  of  figure  12  that  intermittent  breaking  weakens 
the  formation  of  Langmuir  cells  generated  under  uniform  wind  stress.  In  order  to 
expose  this  breaker-vortex-force  interaction  we  show  w -contours  from  four  different 
simulations  with  the  same  Stokes  drift  field  at  two  different  z-levels  (figures  13  and 
14).  These  images  illustrate  an  impact  of  breaking  on  Langmuir-cell  formation  near 
the  water  surface  (0  >  z  >  —2  m)  and  also  surprisingly  deeper  in  the  water  column 
(z  ~  — 13  m).  Overall,  as  the  level  of  intermittent  forcing  increases  (or  equivalently  as 
the  wave  age  decreases)  the  spatial  organization  of  the  Langmuir  pattern  decreases; 
the  lateral  distance  between  downwelling  lines  becomes  wider ;  and  the  lines  shorten  in 
length  in  the  streamwise  direction.  These  trends  are  especially  apparent  if  we  compare 
figure  13(a)  and  figure  13(d),  which  are  simulations  driven  by  uniform  stress  and  by 
large  intermittent  breakers,  respectively.  Also  in  the  simulation  with  the  youngest 
waves  (wave  age  cp/u,a  ~  19)  the  downwelling  is  strongly  focused  at  Y-junctions. 
The  spatial  structure  and  sparsity  of  downwelling  lines  suggest  that  intermittent 
breaking  alters  the  formation  of  Langmuir  cells  compared  to  the  uniform  stress  case. 
The  near-surface  in -patterns  are  consistent  with  the  changes  in  turbulent  variances 
discussed  previously. 

The  flow  visualization  in  figure  14  shows  unexpected  flow  patterns  for  simulations 
with  and  without  breaking  (and  Stokes  drift)  well  below  the  surface  layer,  z  ~  — 13  m. 
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Figure  12.  Snapshot  of  vertical  velocity  uJ  at  z  =  — 1.14  m  for  different  combinations  of 
Stokes  drift  and  wave  breaking  for  fully  developed  waves  cp/u,a  =  30:  (a)  no  wave  effects  and 
uniform  stress,  ( b )  uniform  stress  plus  Stokes  drift,  (c)  breaking  waves  and  no  Stokes  drift, 
and  (d)  breaking  waves  plus  Stokes  drift.  The  grey-scale  bar  shown  at  the  top  of  the  figure  is 
in  units  of  metres  per  second. 


For  flows  driven  by  uniform  stress  or  weak  intermittent  breaking  (figure  14a  and  14 b) 
broad  horizontal  streaks  appear,  and  the  streaks  are  rotated  to  the  right  of  the  surface 
wind  by  Coriolis  effects  (n.b.,  mean  current  profiles  in  figure  5).  These  patterns  are 
signatures  of  depth-filling  Langmuir  circulations  and  are  observed  in  all  LES  of  the 
OBL  driven  by  uniform  stress  and  vortex  force  (e.g.  Skyllingstad  &  Denbo  1995; 
McWilliams  et  al.  1997;  Noh  et  al.  2004).  However,  this  pattern  remarkably  changes 
when  the  surface  forcing  shifts  to  larger  and  more  intermittent  breaking.  The  streaky 
pattern  is  less  organized  as  coherent  round  ‘spots  or  jets’  of  concentrated  negative  w 
appear. 


6.4.  Action  of  vortex  force  on  breaker  vorticity 
The  intriguing  interactions  between  breaking  and  vortex  force  shown  in  figures  13  and 
14,  supported  by  other  visualization,  motivated  a  search  for  the  source  of  downwelling 
jets  in  LES.  Numerous  simulations  with  different  PDFs  of  breaking  show  that  these 
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Figure  13.  Contours  of  downwelling  velocity  w  ^  —0.01  m  s_1  at  z  =  —0.38  m  for  cases  with 
Stokes  drift  and  varying  wave  age:  («)  uniform  stress  no  breaking,  ( b )  fully  developed  waves 
cp/u,a  =  30,  (c)  cp/u,a  =  23,  and  ( d )  cp/u,a  =  19. 


downwelling  jets  are  robust  features.  For  a  given  value  of  Stokes  drift  (or  wind 
speed  and  wave  age)  they  appear  most  often  when  the  forcing  is  dominated  by  large 
intermittent  breakers.  Results  with  a  PDF  that  emphasizes  large-scale  breaking  are 
given  in  Sullivan  et  al.  (2005). 

Our  present  interpretation  for  the  development  of  downwelling  jets  is  based  on  a 
delicate  coupling  of  breaker  vorticity  and  the  so-called  CL2  instability  mechanism 
(Leibovich  1983,  p.  402)  sketched  in  figure  15.  CL2  is  a  potent  pathway  to  generating 
Langmuir  circulations  that  depends  on  positive  feedback  between  the  vortex  force  and 
local  perturbations  of  horizontal  current.  For  example,  Nepf  et  al.  (1995)  speculate 
that  wave  breaking  in  their  channel  flow  may  provide  ‘seed  vorticity’  for  the  initiation 
of  the  CL2  mechanism.  However,  their  results  are  made  ambiguous  by  the  presence 
of  other  transverse  structures  in  their  wavy  flow,  and  the  large  influence  of  the 
bottom  boundary  layer  in  the  shallow  channel  with  surface  wavelengths  comparable 
to,  or  even  larger  than,  the  channel  depth.  Earlier,  Csanady  (1994)  had  speculated 
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Figure  14.  Vertical  velocity  contours  at  z  =  —13.38  m  for  the  same  flows  as  in  figure  13. 
Note  the  appearance  of  the  coherent  round  down  welling  jets  in  (c)  and  (d).  The  grey-scale  bar 
shown  at  the  top  of  the  figure  is  in  units  of  metres  per  second. 


that  surface  stress  anomalies  of  finite  extent,  due  to  either  wind  gusts  or  breaking 
waves,  might  lead  to  a  ‘forced  CL2  mechanism’,  with  Stokes  drift  tilting  the  lines  of 
vertical  vorticity  at  the  edges  of  the  anomaly,  without  the  need  for  the  feedback  of 
the  CL2  mechanism.  The  idealized  sequence  of  CL2  events  (see  figure  15)  assumes  a 
streamwise  current  fluctuation  u'  of  finite  y-extent  and  a  Stokes  drift  uSr  and  requires 
no  coherent  surface-wave  structure.  First,  the  current  anomaly  generates  vertical 
vorticity  of  opposite  signs  that  lead  to  opposing  vortex  forces  uStco^y  that  decay  with 
depth  providing  the  torque  to  generate  counter-rotating  vortices.  Then  the  action  of 
the  vortex  forces  directed  inward  towards  the  centerline  of  the  current  fluctuation 
reinforces  the  perturbation  and  promotes  the  instability.  The  flow  convergence  at  the 
surface  leads  to  downwelling.  In  the  original  description  of  CL2  Leibovich  (1983, 
p.  399)  states  ‘Vorticity  in  the  water  body  may  arise  from  currents  whose  origins 
are  unspecified,  or  ...  by  an  applied  wind  stress  . . .’.  We  find  the  intensity,  scale 
and  proximity  of  the  vorticity  generated  by  breaking  are  important  variables  to 
be  considered  in  CL2.  These  dependences  might  be  anticipated  since  figure  15  is  an 
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Figure  15.  Sketch  illustrating  the  CL2  instability  mechanism  for  the  generation  of  Langmuir 
circulations  at  the  water  surface,  adapted  from  Leibovich  (1983).  The  current  perturbation 
u(y,  z)  generates  vertical  vorticity  coz  that  is  amplified  by  the  vortex  force  V /  =  us,a>y  that  leads 
to  downwelling  at  a  convergence  line.  Breaking  waves  are  a  source  of  the  initial  streamwise 
current  perturbation  and  vertical  vorticity. 

idealization  of  a  single  isolated  perturbation,  and  thus  this  picture  cannot  be  expected 
to  predict  new  dynamics  arising  from  the  coupling  of  nearby  strong  CL2  events. 

To  illustrate  the  importance  of  background  vorticity  on  Langmuir  dynamics  we 
compare  vertical-vorticity  fields  from  simulations  driven  by  uniform  stress  and  by 
intermittent  breaking  in  figure  16.  In  the  former  case  the  vorticity  is  uniformly 
distributed  in  space,  relatively  weak  and  of  small  scale.  In  contrast,  in  the  simulation 
with  breaking  the  vertical  vorticity  is  stronger  in  magnitude  and  contains  a  spectrum 
of  large  and  small  scales.  An  important  point  to  notice  is  that  each  breaker  generates 
both  positive  and  negative  vertical  vorticity  in  close  proximity  that  ideally  interacts 
with  the  Stokes  drift  in  the  vortex  force  to  generate  strong  local  lateral  convergence. 
Inspection  of  the  contours  and  animations  shows  that  as  time  advances  the  initial 
vertical  vorticity  field  from  each  breaker  evolves  into  multiple  pairs  of  plus-minus 
signed  vorticity,  all  of  which  merge  at  a  vigorous  forward  looking  Y-junction;  the 
magnitude  of  downwelling  at  a  Y-junction  with  breaking  exceeds  values  obtained 
under  uniform  stress.  The  dynamics  in  figure  16  gives  rise  to  the  different  patterns  of 
convergence  lines  shown  in  figure  13. 

Previous  studies  have  found  that  breaker  turbulence  can  persist  for  multiple  wave 
periods  (Melville  et  al.  2002;  Sullivan  et  al.  2004)  despite  the  limited  lifespan  of 
active  breaking,  a  time  scale  of  the  order  of  the  wave  period  T  =  2t ic/g.  Thus  it 
is  interesting  to  examine  the  interaction  between  breaking  turbulence  and  the  vortex 
force  for  t  >  T.  The  post-breaking  life  cycle  of  vertical  vorticity,  vertical  velocity 
and  streamwise  current  for  a  typical  breaking  event  extracted  from  a  simulation 
is  displayed  in  figures  17-19.  The  breaker  scales  are  X  ~  20m,  c  ~  5.6 ms-1  and 
T  ~  3.6  s.  The  time  and  space  evolution  of  the  flowfields  is  qualitatively  similar  to 
the  flow  pattern  sketched  in  figure  15  and  illustrates  that  large  breaking  events  seed 
the  CL2  instability.  However,  complex  dynamics  can  occur  at  merger  points.  At  the 


Figure  16.  Snapshot  of  vertical  vorticity  a>-z  at  z  =  — 1.14  m  for  two  simulations  with  Stokes 
drift  driven  by:  (a)  uniform  stress  and  ( b )  breaking  with  wave  age  cp/u,a  =  19.  Note  the  paired 
plus  and  minus  signed  vertical  vorticity  that  occurs  at  the  lateral  (y)  ends  of  each  breaker.  The 
colour  bar  shown  at  the  top  of  the  figure  is  in  units  of  per  second. 


initial  time,  the  residual  of  a  breaking  event  is  a  pair  of  vertically  aligned  vortices 
with  a  y-separation  approximately  equal  to  X  (we  refer  to  these  as  the  primary 
vortices).  Each  of  the  primary  vortices  amplifies  and  migrates  towards  a  common 
centerline  as  a  result  of  the  vortex  force.  With  increasing  t,  current  perturbations 
and  vortices  internal  to  the  primary  pair  develop  as  a  result  of  CL2.  This  appears 
as  an  enhancement  of  the  n-current  at  the  internal  edge  of  each  primary  vortex  (see 
figure  \lb,d).  The  new  four  vortex  system  grows  in  strength,  propagates  forward  and 
heads  to  an  inevitable  collision  at  a  forward  location.  However,  notice  the  sign  of  the 
internal  vortices  is  such  that  they  collaborate  to  generate  strong  backflow  (negative  u ) 
along  the  centerline  of  the  vortex  system  retarding  the  forward  motion  of  the  breaker 
impulse.  The  flow  response  is  vigorous  downward  motion.  At  late  time,  after  about 
80  s  or  22  wave  periods,  the  vortices  all  merge  ending  the  event. 

The  response  of  the  vertical  velocity  field  at  z  =  —1.14  m  shown  in  figure  18  is 
consistent  with  the  active  movement  of  the  vortices.  Immediately  after  the  breaking 
event  w  is  laterally  (y)  distributed  in  front  of  and  behind  the  breaker,  a  response  to 
the  forward  impulse  of  breaking.  As  the  CL2  instability  grows  a  pair  of  downwelling 
lines  develop  just  inside  the  primary  outer  vortices.  The  downwelling  intensifies  until 
at  late  time  the  lines  join  (intense  Y-junctions  are  readily  apparent  in  figure  13d). 
Upstream  of  the  downwelling  merger  the  w-currents  slow  and  change  sign  as  if  in 
response  to  a  stagnation  point. 

Deeper  in  the  water  column,  z  ~  — 13  m,  a  round  downwelling  jet  forms  slightly 
forward  and  to  the  right  of  the  surface  breaking  event.  At  this  depth  the  jet  first 
appears  about  80-100  s  after  the  breaking  event.  Its  appearance  is  consistent  with  a 
travel  time  of  (9(100)  s  and  a  vertical  velocity  w  ~  —0.1  m  s-1.  Thus  we  are  lead  to 
conclude  that  the  downwelling  jets  in  figure  19  (and  those  in  figure  14d)  originate 
at  the  water  surface  and  result  from  interactions  between  CL2  type  instabilities.  We 
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Figure  17.  Temporal  evolution  of  breaker  flowfields  under  the  action  of  the  vortex  force  at 
z  =  —1.14  m  for  the  simulation  with  wave  age  cp/u,a  =  19.  The  resolved  vertical  vorticity  ZH  z 
is  shown  in  the  colour  images  and  the  resolved  streamwise  current  u  as  solid  contour  lines. 
The  colour  bar  at  the  top  of  the  figure  is  in  units  of  per  second  and  the  contour  lines  are 
evenly  spaced  in  the  range  [0,  0.3]  m  s-1.  Relative  to  panel  (a)  the  elapsed  time  in  seconds  is 
[18, (a)],  [37, (Z?)],  [44, (J)],  [55, (e)],  [80,(/)].  Note  the  expanded  v-scale  beginning  with  panel  (c). 


find  that  the  scale  and  intensity  of  the  downwelling  jets  varies  as  expected  with  the 
breaker  phase  speed  c.  With  intermittent  large-scale  breaking  the  downwelling  jets 
penetrate  well  below  the  surface  layer  and  alter  the  structure  and  mixing  of  the  OBL. 

The  flow  patterns  displayed  in  figures  17-19  are  typical  of  strong  breaking  events. 
Breaking  supplies  seed  vertical  vorticity  to  start  the  CL2  instability,  but  the  flow 
can  rapidly  degenerate  into  multiple  pairs  of  vortices  each  acting  under  the  influence 
of  the  vortex  force.  The  forward  propagating  vortex  system  merges  and  produces 
downwelling.  Thus  surface  interactions  between  breaker  vorticity  and  vortex  force 
lead  to  the  formation  of  a  new  coherent  structure  below  the  surface  layer.  When  the 
vertical-vorticity  seeds  of  the  CL2  instability  are  relatively  weak,  streamwise  oriented 
vortices  develop  that  evolve  into  classic  Langmuir  circulations.  However,  if  the  initial 
vorticity  is  strong,  downwelling  jets  also  develop  and  the  surface  Langmuir  pattern 
is  more  intermittent.  Our  results  show  that  turbulence  from  breaking  and  the  CL2 
mechanism  play  a  role  in  an  OBL  with  fully  developed  turbulence,  a  wave  spectrum, 
and  realistic  currents. 
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Figure  18.  Temporal  evolution  of  resolved  vertical  velocity  uJ  and  stream  wise  current  u  for 
the  same  time  periods  and  spatial  locations  as  in  figure  17.  The  colour  bar  at  the  top  of  the 
figure  is  in  units  of  metres  per  second. 


The  vertical  velocity  skewness  (uJ3)/(u>2)3/2  is  a  bulk  statistical  measure  often  used 
to  detect  the  presence  of  coherent  structures  in  a  boundary-layer  flow.  For  example, 
Moeng  &  Rotunno  (1990)  find  that  a  bias  in  vertical  velocity  skewness  indicates  the 
presence  of  thermal  plumes  in  a  convective  flow.  Figure  20  shows  vertical  velocity 
skewness  from  our  simulations  with  different  combinations  of  breaking  and  vortex 
force.  In  all  cases  with  intermittent  breaking  and  vortex  force  the  skewness  is  more 
negative  than  the  comparable  case  driven  by  uniform  stress  and  Stokes  drift.  In 
the  absence  of  the  vortex  force  the  skewness  is  only  slightly  negative.  The  large 
negative  skewness  in  simulation  53  + Sr  reflects  the  presence  of  vigorous  downwelling 
jets  discussed  above.  Notice  in  this  simulation,  the  skewness  is  persistently  negative 
well  below  the  surface  layer  indicating  that  once  the  downwelling  jets  develop  at 
the  surface  they  persist  throughout  the  mixed  layer.  This  is  additional  evidence 
that  strong  intermittent  breaking  and  the  vortex  force  can  act  synergistically  to 
form  a  coherent  structure  roughly  analogous  to  a  thermal  plume.  The  bias  towards 
negative  skewness,  induced  by  surface  waves,  appears  to  be  robust  across  wind  speed 
as  a  simulation  at  Ua  =  30  m  s-1  shows  a  similar  trend,  with  further  discussion 
in  §7. 
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Figure  19.  Time  history  of  vertical  velocity  u;  at  z  =  —13.38  m  for  the  same  flow  conditions 
as  in  figures  17  and  18.  Relative  to  figure  17(a)  the  elapsed  time  in  seconds  is  [44 ,(d)],  [55, (e)], 
[80,(/)],  [90, (g)],  [117,(/i)],  [143,(1)].  The  colour  bar  at  the  top  of  the  figure  is  in  units  of  metres 
per  second.  At  this  depth  the  downwelling  jet  forms  at  nearly  the  same  horizontal  location  as 
the  breaking  event  in  figure  17,  but  begins  to  appear  80-90  s  later. 


6.5.  Scalar  statistics  and  mixed-layer  depth 

The  importance  of  surface  waves  to  OBL  dynamics  and  scalar  transport  is 
traditionally  assumed  to  be  confined  to  a  depth  on  the  order  of  the  wave  height. 
However,  these  LES  results  clearly  indicate  that  surface  waves  alter  mixed-layer 
dynamics  over  a  greater  depth  (9(\h\),  mainly  due  to  the  vortex  force.  Scalar  transport 
is  important  in  the  mixed  layer,  especially  so  for  OBLs  under  high  wind  conditions 
since  entrainment  cooling  can  alter  the  development  of  tropical  cyclones  (e.g.  Emanuel 
1999,  2004).  This  raises  the  important  question  of  whether  surface  waves  play  an 
expanded  role  in  scalar  mixing. 

Figure  21  illustrates  the  influence  of  surface  waves  on  the  bulk  mixed  layer 
temperature  and  the  structure  of  the  thermocline  inversion.  Averaged  over  identical 
periods,  the  mixed-layer  temperature  is  cooler  in  the  presence  of  vortex  force  due 
to  more  efficient  entrainment  at  the  thermocline.  Compared  to  the  cases  driven 
by  uniform  stress  or  breaking  only  the  simulations  with  vortex  force  also  exhibit 


441 


Ocean  boundary  layers  with  vortex  force  and  stochastic  breaking 


(w3)/^)3'2 


Figure  20.  Profiles  of  vertical  velocity  skewness  (w3)  /  (w2)3?2 :  no  wave  effects  dash-dot  line; 
breaking-only  wave  age  cp/u,a  =  30,  solid  line;  Stokes  drift  only,  ►;  Stokes  drift  plus  breaking 
with  wave  age  cp/u,a  =  30,  •;  Stokes  drift  plus  breaking  with  wave  age  cp/u,a  =  23,  o;  Stokes 
drift  plus  breaking  with  wave  age  cp/u,a  =  19,  V;  and  Stokes  drift  plus  breaking  with  wave 
age  cp/u,a  =  30  for  Ua  =  30  m  s_1,  □. 


Figure  21.  Profiles  of  normalized  temperature  for  simulations:  no  wave  effects,  dash-dot  line; 
breaking-only  wave  age  cp/u,a  =  30,  solid  line;  Stokes  drift  only,  ►;  and,  Stokes  drift  plus 
breaking  with  wave  age  cp/u,a  =  19,  V.  The  mixed-layer  depth  and  temperature  are  normalized 
by  the  initial  value  /i;  =  —32.4  m  and  the  reference  temperature  6r  =  283.5  K,  respectively. 
The  time  averaging  is  over  10000  s  starting  at  t  =  40000  s  for  each  simulation.  For  reference, 
the  initial  mixed-layer  sounding  is  shown  as  a  dotted  line.  Notice  how  the  temperature  profiles 
from  all  simulations  relax  back  to  the  initial  sounding  for  z/|A;  <C  —  1. 
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Figure  22.  Profiles  of  vertical  (total)  scalar  flux  normalized  by  the  surface  value  Q,  for 
simulations:  no  wave  effects,  dash-dot  line;  breaking  only  wave  age  cp/u,a  =  30,  solid  line; 
Stokes  drift  only,  ►;  Stokes  drift  plus  breaking  with  wave  age  cp/u,a  =  30,  •;  Stokes  drift 
plus  breaking  with  wave  age  cp/u,a  =  23,  o;  and  Stokes  drift  plus  breaking  with  wave  age 
Cp/u*a  =  19,  V. 


sharper  temperature  gradients  in  the  entrainment  zone.  In  these  simulations  vigorous 
turbulence  generated  by  the  combined  action  of  Langmuir  circulations  and  breaking 
overshoots  the  mean  entrainment  height  leading  to  a  stiffer  inversion;  this  is  analogous 
to  entrainment  dynamics  in  the  daytime  convective  atmospheric  boundary  layer 
(Sullivan  et  al.  1998).  Thus  our  simulations  show  that  surface  waves  can  modulate 
the  structure  of  the  OBL  entrainment  zone.  This  is  further  illustrated  in  figure  22 
where  we  compare  vertical  profiles  of  total  average  scalar  flux  (w9')  for  our  suite 
of  simulations.  Recall,  the  same  small  surface  flux  Q,  is  imposed  in  all  simulations. 
As  expected  for  these  neutrally  stratified  OBLs,  scalar  mixing  is  dominated  by 
entrainment  dynamics.  The  increase  in  entrainment  in  the  simulations  with  the  vortex 
force  shows  that  surface  waves  do  indeed  enhance  mixing  in  the  OBL  compared  to  a 
baseline  simulation  with  no  Stokes  drift.  This  is  largely  a  consequence  of  depth-filling 
Langmuir  cells  that  promote  mixing,  but  also  because  the  Langmuir  cells  alter  the 
structure  of  the  current  profiles  near  the  thermocline  (Skyllingstad  2005),  i.e.  the 
Langmuir  cells  alter  the  current  shear  near  z  ~  h.  The  down  welling  jets  in  simulation 
53  +  St  are  also  observed  to  play  a  critical  role  in  setting  the  scalar  mixing  efficiency 
of  the  OBL;  the  entrainment  flux  is  largest  in  this  simulation.  The  entrainment  of 
cooler  water  from  below  the  thermocline  can  increase  by  almost  a  factor  of  4  in 
the  simulations  with  both  wave  breaking  and  Stokes  drift  compared  to  the  baseline 
simulation  with  no  wave  effects.  At  the  same  time,  surface  waves  enhance  the  scalar 
variance  as  shown  in  figure  23.  The  Langmuir  cells  (and  downwelling  jets)  generated 
by  surface  waves  alter  the  shape  of  the  thermocline  temperature  profile  and  potentially 
are  powerful  enough  to  excite  internal  gravity  waves  (Chini  &  Leibovich  2003).  These 
results  are  not  impacted  by  the  SGS  model  for  scalar  flux.  For  the  simulations  (w6r)sgs 
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Figure  23.  Scalar  variance  profiles  for  simulations  with:  no  wave  effects,  dash-dot  line; 
breaking  only  wave  age  cp/u,a  =  30,  solid  line;  Stokes  drift  only,  ►;  Stokes  drift  plus  breaking 
with  wave  age  cp/u,a  =  30,  •;  Stokes  drift  plus  breaking  with  wave  age  cp/u,a  =  23,  o; 
and  Stokes  drift  plus  breaking  with  wave  age  cp/u,a  =  19,  V.  The  normalizing  parameter 
0,  =  Q,/u,„. 


referenced  to  the  maximum  entrainment  flux  in  figure  22  peaks  at  15%  at  z/\h\  ~  —  1 
and  over  the  bulk  of  the  mixed  layer  —0.9  <  z/\h\  <  0  is  less  than  5%. 


7.  Wind  speed  dependence 

The  preceding  LES  results  for  a  wind  speed  of  15  m  s_1  clearly  demonstrate  that 
OBL  dynamics  and  mixing  depend  on  the  structure  and  state  of  the  surface  wave 
field,  through  the  wave  height  and  breaker  spectra.  These  spectra  vary  with  Ua  (and 
for  a  fixed  wind  also  with  wave  age  cp/u,a),  so  it  is  important  to  consider  how 
the  results  might  change  for  lower  and  higher  winds,  and  in  situations  of  wind-wave 
disequilibrium  as  the  latter  is  the  more  probable  state  of  winds  and  waves.  Future  LES 
will  more  fully  quantify  these  dependences,  but  we  make  some  preliminary  estimates 
assuming  the  empirical  formulas  in  §§  2  and  3  hold,  in  particular  that  the  Froude 
scaling  is  valid  for  different  Ua.  Ideally  for  an  LES  posing  we  need  both  the  \Ja  and 
cpl u*a  dependence  of  the  Stokes  drift,  the  atmospheric  momentum  and  energy  fluxes 
and  the  breaker  spectrum.  We  are  also  interested  in  large  values  of  Ua,  where  imagery 
of  the  sea  surface  shows  impressive  breaking  events,  and  OBL  dynamics  regulating 
the  sea-surface  temperature  are  especially  critical  for  tropical  cyclone  evolution. 

First,  the  relative  importance  of  Stokes  forcing  compared  to  wind  stress  depends 
on  the  turbulent  Langmuir  number  Lat  =  \Ju,aluSt  (McWilliams  et  al.  1997)  with  the 
importance  of  conservative  wave-current  interactions  increasing  with  decreasing  Lat. 
Li,  Garrett  &  Skyllingstad  (2005)  extensively  explore  the  parameter  space  spanned 
by  buoyancy,  shear  and  Stokes  forcing  using  LES,  and  for  typical  oceanic  conditions 
find  the  OBL  is  in  a  regime  dominated  by  Langmuir  turbulence  Lat  ~  0.3,  i.e. 
wave-current  interactions  are  important.  To  test  this  Lat  scaling  in  the  presence  of 
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Figure  24.  The  impact  of  varying  Stokes  drift  on:  (a)  mean  currents,  ( b )  vertical  velocity 
skewness,  (c)  turbulence  variances,  and  ( d )  scalar  flux.  The  solid  line  in  each  figure  is  the 
result  using  the  wave  equilibrium  spectrum  given  by  (2.1)  and  the  dotted  line  is  the  result 
for  a  fetch-limited  wave-height  spectrum  given  by  Komen  et  al.  (1994,  p.  187)  with  wave  age 
cp/u,a  ~  19.  The  magnitude  of  the  Stokes  drift  uSl  for  the  fetch-limited  spectrum  is  reduced 
by  about  a  factor  of  2  compared  to  the  equilibrium  spectrum. 


intermittent  breaking  waves  an  LES  identical  to  63  +  St  is  performed  but  with  a 
reduced  value  of  Stokes  drift  corresponding  to  growing  seas  with  wave  age  cp/u,a  =  19 
(case  B 3  +  FSt  in  table  1).  In  this  LES,  the  Stokes  drift  profile  is  computed  from 
the  empirically  determined  fetch-limited  wave  spectrum  given  by  Komen  et  al.  (1994, 
p.  187).  This  is  a  severe  test  of  the  solution  robustness  as  uSl  is  reduced  by  more  than 
a  factor  of  2  in  the  upper  part  of  the  mixed  layer  when  the  wave  age  varies  from  30 
to  19;  this  change  in  wave  age  implies  La,  increases  by  *Jl.  A  comparison  of  various 
statistical  moments  from  simulations  B 3  +  St  and  B 3  +  FSt  is  given  in  figure  24. 

The  overall  impression  from  these  results  is  that  changing  the  Stokes  forcing  by 
a  factor  of  2  does  not  alter  the  basic  qualitative  behavior  of  the  OBL.  The  shape 
of  the  current  profiles,  enhanced  spanwise  variance  and  large  entrainment  flux  are 
evidence  that  wave-current  interactions  and  breaker  forcing  are  important  to  the 
flow  dynamics  over  a  wide  range  of  sea  states.  Flow  visualization  of  the  ic-current 
(not  presented)  clearly  shows  streamwise  oriented  Langmuir  circulations  and  the 
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formation  of  downwelling  jets  in  simulation  B 3  +  FSt  similar  to  those  in  B 3  +  St. 
The  negative  bias  in  the  vertical  velocity  skewness  in  figure  24  reflects  the  presence 
of  these  coherent  structures.  Although  clearly  not  exhaustive,  our  solutions  appear 
reasonably  robust  to  changes  in  Stokes  forcing  and  further  support  the  Lat  scaling 
proposed  by  McWilliams  et  al.  (1997)  and  verified  by  Li  et  al.  (2005)  in  the  presence 
of  breaking  waves. 

To  investigate  the  wind  speed  dependence  further  we  next  utilize  the  formula  in  §  2, 
u,a  ~  \fOdU a  and  uSr  ~  1  / fp  or  ~  Ua.  Therefore  the  turbulent  Langmuir  number  has 
a  weak  wind  speed  variation  because  of  the  drag  dependence  Lat  ~  Cd25.  For  Ua  = 
15  m  s”1  and  the  Stokes  velocity  at  z  =  1  m,  we  estimate  Lat  =  0.3,  which  is 
a  regime  where  Stokes  forcing  is  important  (see  also  McWilliams  et  al.  1997;  Li 
et  al.  2005).  Doubling  the  wind  speed  to  30  m  s-1  increases  the  drag  coefficient 
Cd  =  2.3  x  10~3  and  therefore  Lat  =  0.33.  At  even  higher  winds,  Ua  >  30  m  s-1, 
Cd  is  probably  wind  speed  independent  or  even  decreasing  (Donelan  et  al.  2004), 
which  again  implies  Lat  ~  0.3.  Based  on  the  above  scaling  we  conclude  for  all  winds 
greater  than  15  m  s_1  the  vortex  force  should  remain  significant  relative  to  the  wind 
stress. 

Estimating  the  high  wind  speed  dependence  of  the  breaker  spectrum  is  far  less 
certain  due  to  the  sparsity  of  field  data.  In  §§2  and  3  we  use  bulk  momentum  and 
energy  conservation  arguments  along  with  data  from  Terray  et  al.  (1996),  Melville  & 
Matusov  (2002),  Melville  et  al.  (2002)  in  the  wind  speed  range  Ua  ~  7—15  m  s”1  to 
build  a  breaker  distribution  with  Ua  and  cp/u,a  dependences.  We  now  assume  these 
scaling  arguments  hold  at  higher  winds.  The  critical  trends  for  the  breaker  distribution 
are  given  in  figures  1-3.  The  peak  phase  speeds  (or  scales)  for  momentum  and  energy 
transfer  (cT,  cE)  from  breaking  move  closer  to  the  peak  in  the  wave-height  spectrum 
with  increasing  wind  speed  and  decreasing  wave  age  (figure  3).  An  estimate  for  the 
peak  scale  for  breaker  momentum  can  be  obtained  by  noticing  that  the  integrand  of 
(3.16)  changes  sign  at  a  crossover  point  c  ~  g,u»a.  Using  this  as  an  estimate  for  cT 
yields 


cT/cp  *  2nvg,u,a/Ua  or  cT/cp  «  g,\/rCd  , 


(7.1) 


which  shows  a  slightly  superlinear  dependence  of  breaker  momentum  on  wind 
speed  (contained  in  cp  and  the  drag  coefficient)  and  linear  dependence  on  wave 
age  (contained  in  the  Terray  parameter). 

Based  on  the  above  arguments  the  general  conclusion  is  that  with  increasing  Ua  the 
vortex  force  remains  roughly  constant  and  the  importance  of  intermittent  breaking 
increases.  Results  from  LES  at  \Ja  =  30  m  s-1,  outlined  in  §5,  support  these  scaling 
arguments.  The  vertical  velocity  skewness  in  figure  20  at  double  the  wind  speed 
has  a  local  minimum  near  the  surface  and  a  broad  region  near  the  thermocline, 
indicating  the  presence  of  coherent  structures  near  z  ~  h.  Flow  visualization  indicates 
downwelling  jets  form  at  a  wave  age  cp/u,a  ~  30  at  this  wind  speed.  (Recall  that 
at  lower  winds  the  downwelling  jets  were  either  quite  weak  or  nonexistent  for  the 
same  value  of  cp/u»a.)  This  is  a  direct  consequence  of  the  shift  to  larger-scale  more- 
intermittent  breaking  shown  in  the  PDF  P{c)  of  figure  1.  Compared  with  a  baseline 
simulation  with  no  wave  effects,  the  entrainment  of  cool  water  at  the  thermocline 
increases  by  at  least  30%  due  to  wave  processes,  and  this  is  likely  to  be  consequential 
for  the  sea-surface  temperature  evolution. 
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8.  Final  remarks 

8.1.  Conclusions 

A  large-eddy  simulation  model  of  the  OBL  is  described  that  accounts  for  surface  wave 
effects  through  non-conservative  breaking  waves  and  phase-averaged  conservative 
wave-current  interactions.  The  equations  for  the  resolved  flow  components  are  the 
Craik-Leibovich  (CL)  theory  with  crucial  wave-current  coupling  through  the  vortex 
force.  These  equations  are  further  augmented  by  a  sum  of  discrete  momentum 
impulses  to  model  breaking.  A  broadband  spectrum  of  waves,  typical  of  measured 
conditions,  is  used  to  compute  the  Stokes  drift  profile  which  appears  in  the  vortex 
force.  The  LES  equations  are  closed  using  an  eddy  viscosity  approach  based  on  the 
subgrid-scale  turbulent  kinetic  energy.  The  prognostic  subgrid-scale  TKE  equation 
includes  Stokes-drift  production  and  breaker  work  terms  that  enhance  subgrid-scale 
energy  near  the  water  surface.  In  this  LES  the  traditional  method  of  forcing  an 
OBL  by  constant  surface  stress  is  replaced  by  a  stochastic  model  that  intermittently 
supplies  momentum  and  energy  to  the  underlying  currents.  The  aggregate  effect  of 
breaking  is  a  linear  superposition  of  discrete  events  of  varying  scale  with  each  event  a 
time-varying  three-dimensional  impulse.  The  breaker  model  obeys  Froude  scaling  and 
is  designed  to  match  laboratory  and  field  observations.  Matching  the  bulk  momentum 
and  energy  fluxes  between  the  atmosphere  and  ocean  determines  constants  in  the  PDF 
of  breaking,  and  the  rate  of  breaking  and  its  variation  with  wind  speed  and  wave  age 
are  critical  components  of  the  model. 

LES  results  at  moderate  wind  speed  Ua  =  15  m  s-1  illustrate  how  surface  waves 
impact  the  dynamics  and  mixing  in  the  OBL  in  important  ways.  Langmuir  circulations 
are  induced  by  the  phase-averaged  wave-current  interactions,  primarily  through  the 
vortex  force.  They  are  depth-filling;  they  increase  the  efficiency  of  vertical  transport; 
and  they  enhance  entrainment  at  the  base  of  the  thermocline.  In  an  OBL  driven  by 
constant  stress  with  a  broad  wave-height  spectrum  only  small  differences  were  found 
from  previous  LES  that  used  a  single  dominant  wave  component.  For  fully  developed 
waves  the  impact  of  intermittent  breaking  waves  is  mainly  confined  to  the  region 
close  to  the  water  surface;  they  energize  the  surface  layer,  destroy  traditional  wall 
layer  scaling  and  enhance  the  turbulent  dissipation. 

Stochastically  distributed  breaking  waves  also  disrupt  the  surface  downwelling 
patterns,  the  characteristic  signatures  of  Langmuir  circulations.  As  the  breaking 
becomes  more  intermittent  and  focused  towards  larger  scales,  the  downwelling  lines 
shorten  in  length  and  become  more  widely  separated  in  the  lateral  direction.  When 
downwelling  lines  join,  they  are  sites  of  concentrated  negative  vertical  velocity.  For 
a  given  wind  speed,  fully  developed  equilibrium  waves  break  often  at  small  scales, 
while  younger  seas  break  more  intermittently  at  larger  scales  closer  to  the  peak 
in  the  wave-height  spectrum.  This  wave-age  dependence  has  consequences  for  the 
interaction  between  breaking  and  Langmuir  circulations.  Although  the  vortex  force 
and  breaking  are  added  linearly  to  the  equations,  they  can  couple  in  a  surprising 
way  to  generate  a  new  coherent  structure,  a  downwelling  jet,  which  weakens  but 
coexists  with  traditional  Langmuir  cells  (i.e.  in  the  absence  of  breakers).  The  origins 
of  downwelling  jets  are  surface  breaking  sites,  and  their  formation  is  closely  linked 
to  the  strength  and  scale  content  of  the  near-surface  vertical  vorticity  field  and  the 
CL2  instability  mechanism.  Streamwise  breaking  and  its  associated  strong  vertical 
vorticity  interacts  with  the  vortex  force  over  multiple  wave  periods.  Breaking  seeds 
the  CL2  instability  and  generates  strong  lateral  convergence  at  a  location  forward  of 
the  initial  breaking  site.  The  onset  of  downwelling  jets  coincides  with  a  shift  toward 
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larger-scale  breaking,  i.e.  as  the  wind  speed  increases  and/or  the  wave  age  decreases. 
The  jets  are  found  to  extend  throughout  the  mixed  layer  and  can  significantly  alter 
the  vertical  velocity  skewness.  Wave  processes  are  globally  important  as  Langmuir 
circulations  and  downwelling  jets  enhance  the  entrainment  of  cool  water  at  the  OBL 
therm ocline.  In  summary,  OBLs  are  importantly  different  from  their  wall-bounded 
shear-layer  counterparts  because  of  surface  waves.  These  effects  need  to  be  accounted 
for  in  bulk  parameterizations  of  the  OBL. 


8.2.  Discussion 

There  are  several  aspects  of  the  LES  predictions  that  invite  further  numerical 
evaluation  and  experimental  testing.  First,  we  adopt  the  asymptotic  CL  model  of 
wave-current  interactions  with  no  feedback.  Fully  coupled  current  and  wave  fields 
might  introduce  transverse  inhomogeneity  that  can  potentially  alter  the  Stokes  drift 
and  breaking  distributions  which  are  used  in  the  present  model.  Work  by  Zhou 
(1999)  and  Kawamura  (2000)  indicates  that  CL  is  a  valid  approximation  but  further 
exploration  is  warranted.  Next,  the  LES  wave  forcing  and  current  response  predictions 
need  to  be  tested  in  field  conditions  for  a  variety  of  wind  speeds  and  wave  states. 
This  has  the  advantage  of  validating  the  LES,  but  it  can  also  help  clarify  the 
PDF  of  breaking  and  its  dependence  on  environmental  conditions.  The  LES  results 
suggest  wave  processes  are  important  at  high  winds,  especially  so  for  tropical  cyclone 
evolution  through  enhanced  entrainment  of  cool  water  at  the  thermocline.  This  needs 
to  be  checked  with  new  instrumentation  and  focused  field  campaigns  and  additional 
explorations  of  the  LES  parameter  space.  The  present  LES  results  are  an  impetus 
to  incorporate  more  fully  wave  influences  in  bulk  models  of  the  OBL  such  as  the 
K-profile  parameterization  (Large  et  al.  1995;  McWilliams  &  Sullivan  2000;  Smyth 
et  al.  2002;  McWilliams  &  Huckle  2006). 

Perhaps  the  greatest  need  is  to  improve  the  description  of  the  breaking  statistics  as 
a  function  of  wind  speed,  wave  age  and  wind-wave  misalignment.  The  statistical 
description  of  breaking  used  here  is  speculative,  and  may  be  criticized  for  its 
incompleteness,  but  we  believe  that  the  significance  of  the  new  dynamics  that  is 
introduced  to  OBL  modeling  by  the  explicit  inclusion  of  both  breaking  and  Langmuir 
turbulence  justifies  the  sometimes  crude  assumptions  made  here.  With  very  recent 
developments  in  both  field  measurements  and  modeling  of  surface  waves  and  wave 
breaking,  we  expect  that  one  of  the  next  developments  will  be  to  use  the  output 
of  wind-wave  numerical  models  to  describe  the  breaking  statistics  for  the  OBL,  i.e. 
the  coupling  of  OBL  and  wind-wave  numerical  models.  We  believe  that  the  present 
results  will  provide  guidance  on  addressing  these  developments. 
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Appendix.  LES  code  details 

A.l.  Numerical  method  with  breakers 

The  base  set  of  LES  equations  for  the  resolved  fields,  subgrid-scale  energy  and  pressure 
described  in  §4  are  solved  using  a  mixed  finite-difference  pseudospectral  scheme 
with  a  third-order  Runge-Kutta  time  stepping  with  a  fixed  Courant-Friedrich-Lewy 
(CFL)  number.  Solution  variables  are  staggered  in  the  vertical  direction:  (u,v,9,p) 
are  located  at  cell  centers  while  (w,  e)  are  positioned  at  cell  faces.  Periodic  lateral 
boundary  conditions  are  assumed  while  a  radiation  boundary  condition  is  imposed  at 
the  lower  boundary  ( z  — *•  —  oo)  of  the  computational  domain  (Klemp  &  Durran  1983). 
The  code  is  parallelized  with  the  Message  Passing  Interface  using  vertical  domain 
decomposition.  The  elliptic  Poisson  equation  for  pressure  is  solved  with  a  custom 
built  direct  solver  utilizing  a  parallel  matrix  transpose. 

The  simulation  code  is  extensively  modified  to  account  for  stochastic  wave  forcing. 
Observations  of  breaking  waves  in  the  equilibrium  wind-wave  regime  report  that  the 
phase  speed  of  detectable  breakers  varies  from  c  «  1  to  c  «  14  m  s_1  for  wind  speeds 
7<  Ua  <  15  m  s_1  (Melville  &  Matusov  2002).  Based  on  the  linear  dispersion  relation 
the  breaker  length  scale  k  then  ranges  from  0.6  m  to  more  than  100  m,  while  its  period 
T  ranges  from  0.6  s  to  more  than  9  s.  Capturing  the  breaker  forcing  adequately  over 
this  wide  scale  range  requires  careful  evaluation  of  the  impulse  A,  in  our  integration 
scheme.  From  preliminary  solutions  we  found  that  simply  evaluating  A,(x,  t)  and  then 
inserting  it  into  the  right-hand  side  of  our  standard  algorithm  would  completely  miss 
or  generate  only  a  partial  impulse.  Hence,  as  the  current  fields  moved  forward  in 
time  they  would  either  remain  laminar  or  never  become  fully  turbulent.  The  method 
successfully  employed  advances  the  current  fields  from  stage  n  to  stage  n  +  1  in  the 
Runge-Kutta  scheme  using  the  rule, 

m'!+1(x)  =  u”{x)  +  Atrj  ^r,(x)  -  ^-(x)  +  A*(x,  t^j  4 - ,  (Al) 

where  At  is  the  time  step,  r,  denotes  all  other  terms  on  the  right-hand-side,  of  (4.1), 
?7  is  a  weight  associated  with  the  time  stepping,  and  the  dots  indicate  prior  substeps. 
In  (Al),  Aj  is  the  average  breaker  impulse  defined  as  a  time  and  space  integral  of 
the  breaker  function  (3.11)  averaged  to  the  LES  gridpoint  (x,  t).  For  example,  in  the 
x-component  of  the  current  equations, 

t+At  z+Az/2  y+Ay/2  x+Ax/2 

/  /  /  s/(a,  ft,  S,  y)  cos®  dx  dy  dc  df ,  (A2) 

t  z—Az/2  y—Ay/2  x—Ax/2 

where  the  averaging  is  over  the  small  volume  AV  =  Ax  Ay  Az  and  time  increment 
At,  and  0  is  the  breakers  horizontal  orientation  in  the  LES  grid.  The  integrand 
contains  complex  dependences  on  (x,  t )  because  of  the  breaker  functional  forms 
(3.11)  and  (3.12)  and  thus  the  integrals  in  (A2)  are  evaluated  using  a  trapezoidal 
rule  on  a  fine  mesh  between  LES  gridpoints.  The  implementation  evaluates  (A2)  at 
every  substep  at  each  LES  gridpoint  (x,  t)  for  all  active  breaking  events.  Tests  of  the 
numerical  implementation  show  that  it  conserves  the  total  momentum  from  a  breaker 
irrespective  of  its  time  and  space  scale  and  hence  even  small  scale  breakers  are  ‘felt’ 
in  the  LES  grid. 

In  our  simulations  the  maximum  allowable  timestep  At  is  picked  based  on  the 
grid  resolution,  the  maximum  velocity  in  the  computational  domain  and  a  fixed 
CFL  number.  This  method  is  adaptive  at  each  timestep  and  is  a  robust  scheme 
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for  determining  an  optimum  timestep  for  a  variety  of  flows.  Simulation  diagnostics 
show  breaking  impacts  the  numerical  stability  of  the  code.  Large  breakers  boost  the 
magnitude  of  the  work  term  W  in  the  SGS  e-equation  (4.1),  which  in  turn  locally 
enhances  the  SGS  eddy  viscosity  for  momentum  and  scalars  (ly,  vh).  To  maintain 
numerical  stability  a  minimum  timestep  is  chosen  to  satisfy  a  CFL  constraint  or  a 
viscous  stability  criterion  based  on  the  subgrid-scale  eddy  viscosities  and  grid  spacing, 
e.g.  At  <  Az2 /2vh  (Ferziger  &  Peric  2002),  whichever  is  more  stringent. 

A.2.  Algorithm  outline 

A  synopsis  of  the  algorithmic  components  of  our  OBL  LES  is  the  following: 

(1)  Select  a  wind  speed  Ua. 

(a)  Compute  the  available  atmospheric  momentum  and  energy  using  bulk 
formulas  (2.5)  and  (2.7). 

( b )  Integrate  the  wave  spectrum  (2.1)  to  generate  a  vertical  profile  of  the  Stokes 
drift  (2.4). 

(c)  Compute  the  breaker  PDF  and  rate  of  breaker  generation  N  using  (3.16)  and 
(3.17),  respectively. 

(d)  Based  on  the  LES  grid  resolution,  partition  the  atmospheric  momentum  and 
energy  fluxes  into  resolved  and  subgrid  breakers  as  in  (4.6)  and  (4.7). 

(2)  Integrate  the  LES  equations  (4.1)  forward  in  time.  Breakers  obey  the  following 
rules  at  each  full  Runge-Kutta  time  step: 

(a)  Draw  breaker  speeds  from  the  PDF  and  use  the  linear  dispersion  relation  to 
establish  other  event  properties. 

( b )  Candidate  breakers  are  sorted  by  size  so  that  the  largest  breakers  are  added 
to  the  water  surface  first.  Small  breakers  fill  in  the  remaining  empty  voids. 

(c)  Breaking  events  are  initiated  at  non-overlapping  surface  sites  drawn  from  a 
uniform  random  distribution. 

(d)  Generate  the  four-dimensional  breaker  impulse  functions  for  each  c  and 
average  the  breaker  impulses  to  the  LES  grid.  Add  breaker  acceleration  and 
work  terms  to  the  LES  governing  equations. 

(e)  Maintain  and  update  the  link  list  of  breaking  events  allowing  old  breakers  to 
die,  existing  breakers  to  evolve,  and  new  breakers  to  be  born. 
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ABSTRACT 

Winds  and  waves  in  marine  boundary  layers  are  often  in  an  unsettled  state  when  fast-running  swell 
generated  by  distant  storms  propagates  into  local  regions  and  modifies  the  overlying  turbulent  fields.  A 
large-eddy  simulation  (LES)  model  with  the  capability  to  resolve  a  moving  sinusoidal  wave  at  its  lower 
boundary  is  developed  to  investigate  this  low-wind/fast-wave  regime.  It  is  used  to  simulate  idealized  situ¬ 
ations  with  wind  following  and  opposing  fast-propagating  waves  (swell),  and  stationary  bumps.  LES  predicts 
momentum  transfer  from  the  ocean  to  the  atmosphere  for  wind  following  swell,  and  this  can  greatly  modify 
the  turbulence  production  mechanism  in  the  marine  surface  layer.  In  certain  circumstances  the  generation 
of  a  low-level  jet  reduces  the  mean  shear  between  the  surface  layer  and  the  PBL  top,  resulting  in  a  near 
collapse  of  turbulence  in  the  PBL.  When  light  winds  oppose  the  propagating  swell,  turbulence  levels 
increase  over  the  depth  of  the  boundary  layer  and  the  surface  drag  increases  by  a  factor  of  4  compared  to 
a  flat  surface.  The  mean  wind  profile,  turbulence  variances,  and  vertical  momentum  flux  are  then  dependent 
on  the  state  of  the  wave  field.  The  LES  results  are  compared  with  measurements  from  the  Coupled 
Boundary  Layers  Air-Sea  Transfer  (CBLAST)  field  campaign.  A  quadrant  analysis  of  the  momentum  flux 
from  CBLAST  verifies  a  wave  age  dependence  predicted  by  the  LES  solutions.  The  measured  bulk  drag 
coefficient  CD  then  depends  on  wind  speed  and  wave  state.  In  situations  with  light  wind  following  swell,  CD 
is  approximately  50%  lower  than  values  obtained  from  standard  bulk  parameterizations  that  have  no  sea 
state  dependence.  In  extreme  cases  with  light  wind  and  persistent  swell,  CD  <  0. 


1.  Introduction 

An  outstanding  question  in  wind-wave  interaction 
studies  is  the  effect  of  fast-running  waves  or  swell  on 
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the  winds  and  turbulence  in  the  atmospheric  planetary 
boundary  layer  (PBL).  Swell-dominated  wave  fields  oc¬ 
cur  after  the  passage  of  storm  fronts,  propagate  long 
distances  without  significant  dissipation  (e.g.,  see  esti¬ 
mates  in  Cohen  and  Belcher  1999),  and  often  dominate 
the  local  wave-height  variance.  In  this  situation  it  is 
difficult  to  measure  and  isolate  the  contributions  of  lo¬ 
cally  generated  wind  waves  to  the  surface  roughness 
and  stress.  The  impact  of  swell  on  surface  drag  param¬ 
eterizations  is  a  closely  related  question.  Donelan  and 
Pierson  (1987)  and  Donelan  et  al.  (1997)  suggest  that 
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swell  influences  are  strong  and  that  wind-swell  align¬ 
ment  is  an  important  factor  for  the  measured  drag  co¬ 
efficients  (e.g.,  they  report  that  the  drag  increases  by  a 
factor  of  3  for  swell  opposing  the  wind).  Thus,  the  sur¬ 
face  stress  likely  depends  on  wind  speed,  wave  age,  and 
swell  amplitude  and  direction. 

In  a  pioneering  study,  Harris  (1966)  first  reported  the 
formation  of  a  wave-driven  wind  in  a  laboratory  wave 
tank.  Since  then  a  growing  body  of  experimental  evi¬ 
dence  has  documented  unique  marine  surface-layer  dy¬ 
namics  in  the  presence  of  swell:  development  of  low- 
level  jets  (Holland  et  al.  1981;  Miller  1999),  positive 
upward  momentum  flux  (Grachev  and  Fairall  2001; 
Smedman  et  al.  1994,  1999),  mean  velocity  profiles  de¬ 
creasing  with  increasing  height  (e.g.,  Rutgersson  et  al. 
2001),  reduced  turbulence  levels  (Drennan  et  al.  1999), 
and  misalignment  between  surface  stresses  and  mean 
winds  (Grachev  et  al.  2003).  Time  series  of  surface- 
layer  winds  collected  from  the  Research  Platform  Float¬ 
ing  Instrument  Platform  (R/P  FLIP),  reported  by  Miller 
(1999,  p.  122),  clearly  show  the  hourly  transition  from  a 
logarithmic  to  nearly  uniform,  near-surface  wind  profile 
after  a  storm  passage;  coincident  with  the  wind-profile 
change  is  a  rapid  reduction  in  the  turbulent  momentum 
flux.  These  features  appear  to  be  signatures  of  a  wave- 
driven  surface  layer  and  invalidate  the  use  of  Monin- 
Obukhov  similarity  theory  that  often  is  used  to  predict 
air-sea  fluxes  (e.g.,  Rutgersson  et  al.  2001).  The  overall 
impact  of  swell  throughout  the  PBL  contradicts  the  com¬ 
mon  view  that  the  depth  of  the  wave  boundary  layer 
(WBL;  i.e.,  the  region  directly  impacted  by  waves)  is  quite 
shallow,  z  <  3  m  (e.g.,  Makin  and  Mastenbroek  1996). 

The  goals  of  this  study  are  to  develop  and  use  tur¬ 
bulence-resolving  large-eddy  simulations  (LES)  to  im¬ 
prove  the  understanding  of  the  interactions  between 
atmospheric  turbulence  and  surface  waves,  and  to  aid 
in  the  interpretation  of  observations  from  the  Coupled 
Boundary  Layers  Air-Sea  Transfer  (CBLAST  low 
wind)  field  campaign  (Edson  et  al.  2007).  Our  use  of 
LES  to  examine  the  low-wind/fast-wave  regime  in  an 
atmospheric  PBL  is  new  (see  Sullivan  et  al.  2004, 
2006b)  but  we  note  that  Reynolds-average  closure 
models  have  previously  been  used  to  study  some  of  the 
impacts  of  fast-moving  waves  on  marine  surface  layers 
(e.g.,  Gent  and  Taylor  1976;  Gent  1977;  Li  1995; 
Kudryavtsev  and  Makin  2004).  This  LES  study  extends 
our  previous  direct  numerical  simulations  over  a  wavy 
lower  boundary  (Sullivan  et  al.  2000;  Sullivan  and 
McWilliams  2002). 

2.  CBLAST  field  campaign 

Motivation  for  the  present  investigation  stems  from 
observations  collected  in  CBLAST  and  similar  field 


studies;  an  overview  of  the  CBLAST  goals,  measuring 
platforms,  datasets,  and  preliminary  analysis  is  given  by 
Edson  et  al.  (2007).  CBLAST  was  a  major  field  cam¬ 
paign  designed  to  investigate  boundary  layer  processes 
that  couple  the  atmosphere,  wave  field,  and  ocean  un¬ 
der  a  variety  of  low  to  moderate  wind  conditions.  The 
site  for  CBLAST  is  the  Atlantic  Ocean  south  of  Mar¬ 
tha’s  Vineyard,  Massachusetts,  and  intensive  observa¬ 
tion  periods  occurred  in  the  summers  of  2001  and  2003. 
The  output  from  this  field  campaign  is  a  large  observa¬ 
tional  database  gathered  over  multiple  months  using  a 
variety  of  sensors  and  measuring  platforms  on  both 
sides  of  the  air-sea  interface.  One  of  the  unique  mea¬ 
suring  platforms  specifically  developed  for  CBLAST  is 
the  Air-Sea  Interaction  Tower  (ASIT)  shown  in  Fig.  1. 
The  ASIT  is  a  low-profile,  fixed  structure  that  mini¬ 
mizes  flow  distortion  and  removes  the  need  for  motion 
correction.  It  is  exposed  to  effectively  infinite  fetch  for 
south-southwesterly  wind  directions.  Atmospheric  sen¬ 
sors  at  fixed  heights  and  on  a  vertical  profiler  provided 
direct  turbulence  flux  measurements,  wind  profiles,  and 
surface  wave  information.  For  our  purposes  we  use  a 
small  subset  of  the  CBLAST  data,  focusing  on  the  ob¬ 
servations  of  the  surface-layer  winds  and  wave  fields 
gathered  from  the  ASIT. 

Wind-wave  equilibrium  is  the  asymptotic  state  of 
aligned  winds  and  waves  where  the  wave  spectrum  is 
fully  developed  and  the  peak  frequency  and  shape  of 
the  wave  spectrum  are  only  changing  slowly  with  time; 
it  occurs  most  often  at  moderate  to  high  winds  Ua  >  10 
m  s-1.  A  bulk  measure  of  wind-wave  equilibrium  is 
when  the  ratio  of  the  peak  frequency  (or  peak  phase 
speed  Cp )  of  the  wave-height  spectrum  to  a  reference 
atmospheric  wind  Ua  attains  a  limiting  value,  CpIUa  = 
1.2  (e.g.,  Alves  et  al.  2003).  The  CBLAST  observations 
strongly  emphasize  the  nonequilibrium  and  variable  na¬ 
ture  of  winds  and  waves  at  low  winds.  Edson  et  al. 
(2007)  finds  that  the  histogram  of  surface  wind  speed 
collected  over  many  months  exhibits  a  maximum  be¬ 
tween  4  and  6  m  s-1,  with  a  few  excursions  up  to  10 
m  s-1,  and  with  a  highly  preferred  wind  direction  from 
the  southwest,  about  225°  from  north  [see  Figs.  4  and  5 
from  Edson  et  al.  (2007)].  Figure  2  shows  that  the  rela¬ 
tive  orientation  between  U„  and  C„,  that  is,  the  wind- 
wave  angle  f>,  is  nearly  randomly  distributed.  The  his¬ 
togram  of  the  wind-wave  angle  has  a  modest  peak  in 
the  range  of  aligned  winds  and  waves  —30°  <  <fi  <  30° 
but  crossing  winds  and  waves,  <f>  =  ±90°,  and  even 
opposing,  <f>  =  ±180°,  are  equally  probable.  Given  the 
observed  preferred  wind  direction  reported  by  Edson 
et  al.  (2007),  the  variability  in  <f>  must  largely  be  a  con¬ 
sequence  of  nonlocal  wave  components,  that  is,  a  result 
of  remotely  generated  swell  propagating  into  the  ob- 
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Fig.  1.  The  Air-Sea  Interaction  Tower  with  twin  masts  deployed  during  the  OBLAST  field 
campaign.  Sonic  anemometers  mounted  on  the  left  (forward)  mast  translate  vertically  to 
obtain  fine  spatial  resolution  of  the  mean  velocity  and  scalar  profiles,  while  fixed  sonic  an¬ 
emometers  attached  to  the  right  (rearward)  mast  are  used  to  measure  vertical  (turbulence) 
fluxes  of  momentum  and  scalars. 


servation  region.  This  conclusion  is  quantified  by  the 
histogram  of  the  wave  age  parameter  CpIUacos<j) 
shown  in  Fig.  2.1 

The  CBLAST  surface  wind  and  wave  fields  are  found 
to  be  dominated  by  relatively  fast-moving  waves  (or  old 
seas)  I Cp/[/acos<(>l  >  1.2.  The  likelihood  of  wave  age 
lying  outside  the  interval  [0, 1.2]  is  about  75%.  We  note 
that  some  of  the  large  wave  age  values  are  a  conse¬ 
quence  of  crossing  winds  and  waves  when  <f>  =  ±tt/2. 
The  overall  conclusion  remains,  however,  that  when  the 
winds  are  light  the  wave  field  is  most  often  in  disequi¬ 
librium  with  the  local  winds.  Churchill  et  al.  (2006), 
using  the  method  described  by  Hanson  and  Phillips 
(2001),  provides  a  detailed  description  of  the  complex 
CBLAST  wave  fields. 

Surface  waves  are  the  primary  source  of  roughness 


1  Several  definitions  of  wave  age  are  used  in  the  literature  (e.g., 
see  Komen  et  al.  1994;  Alves  et  al.  2003;  Plant  1982).  The  present 
definition  is  adopted  since  it  accounts  for  the  directional  align¬ 
ment  between  winds  and  waves.  In  particular,  the  definition  cap¬ 
tures  the  variability  in  wave  age  and  the  occurrence  of  counterseas 
seen  in  our  data. 


Fig.  2.  Frequency  histogram  of  (top)  wind-wave  angle  <f>  and 
(bottom)  wave  age  CpIUacos(l>  during  CBLAST  for  all  wind-wave 
conditions.  In  the  bottom  panel  the  solid  line  is  the  cumulative 
probability  sum  1  -  So  P(x')  dx' ,  where  p(x)  is  the  probability 
density  function. 


Cumulative  Probability 
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Fig.  3.  Drag  coefficients  obtained  from  three  measurement  levels  (squares,  diamonds, 
circles)  =  (4.0,  6.0,  10.0)  m  during  CBLAST  (Edson  et  al.  2006);  CD  is  referenced  to  a  10-m 
height  and  neutral  conditions.  The  TOGA  COARE  3.0  parameterization  is  indicated  by  the 
dashed  line.  Note  the  negative  values  of  CD  and  increase  in  variability  at  low  winds. 


for  the  PBL,  and  it  is  a  long-standing  goal  of  marine 
surface-layer  research  to  quantify  the  drag  of  the  sea 
surface  as  a  function  of  wind  speed.  Figure  3  shows  the 
variation  of  the  neutral  drag  coefficient  CD  as  function 
of  the  reference  atmospheric  wind  speed  Ua  at  a  height 
z  =  10  m  (Edson  et  al.  2007).  The  results  are  from  three 
different  measurement  heights  along  the  ASIT  gath¬ 
ered  over  the  entire  CBLAST  observation  period.  The 
average  CD  varies  linearly  with  wind  speed  in  close 
agreement  with  the  Tropical  Ocean  and  Global  Atmo¬ 
sphere  Coupled  Ocean-Atmosphere  Response  Experi¬ 
ment  (TOGA  COARE)  3.0  parameterization  (Fairall 
et  al.  2003).  However,  notice  the  largest  scatter  in  CD  is 
at  low  wind  speed  and  that  over  certain  periods,  CD  < 
0.  We  hypothesize  some  of  this  scatter  is  attributed  to 
the  nonequilibrium  state  of  winds  and  waves  at  low 
winds  shown  in  Fig.  2. 

3.  PBL  wavy-surface  problem  formulation 

The  CBLAST  observational  results  illustrate  that  in 
low  to  moderate  winds  the  most  common  state  of  the 
marine  boundary  layer  is  disequilibrium  between  the 
surface-layer  winds  and  the  underlying  wave  field;  usu¬ 
ally  the  wave  field  is  propagating  faster  than  and  at  an 
angle  to  the  mean  surface  wind.  Our  study  is  intended 
to  elucidate  the  impact  of  nonequilibrium  waves  on  tur¬ 
bulence  in  the  low-wind  atmospheric  surface  layer  and 
more  generally  the  PBL.  To  investigate  this  low-wind/ 
fast-wave  regime  a  large-eddy  simulation  model  of  the 
PBL  with  the  ability  to  resolve  moving  sinusoidal 
modes  at  its  lower  boundary  was  developed.  A  descrip¬ 


tion  of  the  LES  model  including  the  governing  equa¬ 
tions,  grid  generation,  and  numerical  method  is  pro¬ 
vided  in  the  appendix.  This  LES  model  is  idealized,  as 
a  complete  simulation  of  turbulent  winds  and  a  fully 
interacting  wave  field  at  all  scales  of  interest  is  not  com¬ 
putationally  feasible.  The  design  of  our  PBL  wavy- 
surface  numerical  experiments  instead  focuses  on  for¬ 
mulating  process  studies  to  expose  impacts  of  fast- 
moving  waves  on  surface-layer  winds. 

Here  we  emphasize  neutrally  stratified  (zero  surface 
heat  flux)  PBLs  where  the  wave  propagation  speed  c  is 
large  compared  to  the  surface  wind  and  also  when  the 
waves  are  stationary,  c  =  0.  Three  regimes  are  consid¬ 
ered,  namely,  wind  following  waves,  wind  opposing 
waves,  and  stationary  bumps  (or  small  hills).  These 
cases  serve  to  illustrate  the  importance  of  wave  phase 
speed  relative  to  wind  speed  (i.e.,  wave  age)  and  the 
orientation  of  winds  and  waves. 

In  our  LES  experiments,  the  imposed  surface  wave  is 
two-dimensional  (i.e.,  has  only  x-z  variations)  with 
wavelength  A  =  100  m,  amplitude  a  =  1.6  m,  low  wave 
slope  lairlX  —  0.1,  and  propagates  in  either  the  positive 
or  negative  x  direction.  Based  on  the  linear  dispersion 
relationship  c2  =  gA/2ir,  the  moving  wave  has  phase 
speed  c  =  12.5  m  s-1.  We  set  the  surface  roughness 
za  =  2  X  10~4  m,  a  typical  value  for  a  low-wind  marine 
boundary  layer  (Donelan  1998;  Fig.  2).  Essentially,  the 
z0  parameterization  accounts  for  the  drag  of  unresolved 
small-scale  waves  riding  on  the  larger-scale  resolved 
swell.  The  geostrophic  winds  are  [( Ug ,  Vg)  =  (5,  0)] 
ms-1  and  the  surface  heat  flux  Q.,  =  0.  Relatively  shallow 
PBLs  are  simulated  with  an  initial  depth  zt  =  400  m. 
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Table  1.  Simulation  properties. 


Run  name 

Ug  (m  s  r) 

£?*  (K  ms  ') 

</>  (°) 

CD  (X103) 

Comments 

ZN 

5 

0.0 

6.8 

1.2 

Flat  z0  surface,  neutral  flow 

BN 

5 

0.0 

19.0 

3.7 

Stationary  bumps,  neutral  flow 

FN 

5 

0.0 

-2.3 

-0.12 

Wind  following  waves,  neutral  flow 

FC 

5 

0.01 

-2.1 

-0.1 

Wind  following  waves,  slight  convection 

FN2 

2 

0.0 

-6.1 

-2.2 

Wind  following  waves,  neutral  flow 

ON 

5 

0.0 

172.3 

6.6 

Wind  opposing  waves,  neutral  flow 

The  initial  temperature  sounding  ddldz  =  0  for  0  <  z  < 
Zi,  with  a  strong  stable  inversion  where  ddldz  =  0.01 
Km1  above  z,-.  In  all  runs  the  Coriolis  parameter  /  = 
10~4  s_1.  Thus  our  numerical  experiments  are  an  ide¬ 
alization  of  a  near-neutral  PBL  above  remotely  gener¬ 
ated  swell  with  light  winds  as  the  LES  wave  age  param¬ 
eter  clUg  >  2. 

To  explore  the  sensitivity  of  the  LES  solutions  we 
also  compare  flows  where  we  add  a  small  amount  of 
surface  heating  <2*  =  0.01  K  m  s_1,  set  the  wave  am¬ 
plitude  a  —  0  to  generate  a  traditional  flat  z0-surface 
PBL,  and  reduce  the  geostrophic  wind  ( Ug,  Vg)  =  (2,  0) 
ms-1.  A  summary  of  the  LES  experiments  discussed 
here  is  provided  in  Tabic  1.  For  easy  identification  the 
table  includes  an  abbreviated  run  name  with  comments 
describing  the  bulk  condition  of  the  simulation,  for  ex¬ 
ample,  case  FN  is  a  simulation  with  wind  following 
waves  and  neutral  stratification  (Q,:  =  0).  Also  in  this 
table,  <j)  is  the  angle  between  the  wave  propagation 
direction  and  the  surface  wind,  and  CD  is  the  drag  co¬ 
efficient  deduced  from  the  LES  data  at  z  =  15  m.  Other 
parameters  in  Table  1  are  as  defined  above. 

For  the  suite  of  LES  experiments,  the  computational 
domain  (1200,  1200,  800)  m  is  discretized  using  (250, 
250,  96)  grid  points  with  the  horizontal  resolution 
nearly  constant  in  physical  space  (Ax,  Ay)  =  (4.8,  4.8) 


m.  As  a  result  of  the  horizontal  mesh,  the  waveforms 
imposed  at  the  lower  boundary  are  well  resolved,  ap¬ 
proximately  25  grid  points  per  wave.  In  x-z  planes  a 
conformal  surface-fitted  mesh  is  constructed  between 
the  surface  wave  zb  =  h(x,  t)  =  a  cos  [ k(x  —  cf)]  and  the 
horizontally  flat  top  of  the  computational  domain  z  = 
zL  (see  the  appendix  and  Fig.  4).  The  vertical  spacing  is 
varied  in  order  to  capture  the  different  scales  of  motion 
in  the  PBL;  tight  spacing  Az  «  1  m  is  used  near  the 
surface,  and  the  spacing  expands  smoothly  with 
Az  ^  14  m  at  z„  to  Az  =  28  m  well  above  the  PBL  at  the 
top  of  the  computational  domain.  Approximately  75 
grid  levels  arc  located  between  the  surface  and  the  PBL 
inversion.  Within  the  PBL  the  grid  aspect  ratio  Ax/Az 
varies  from  about  4.8  at  the  surface  to  about  0.34  at  the 
PBL  inversion.  At  the  surface  the  mesh  aspect  ratio  is 
just  within  the  acceptable  limits  for  anisotropic  LES 
grids  (Scotti  et  al.  1993).  The  LES  are  computationally 
intensive.  The  time  step  is  limited  by  the  Courant- 
Fredrichs-Lewy  (CFL)  condition  based  on  the  speed  of 
the  wave  c  and  the  fine  horizontal  spacing  Ax;  typical 
time  steps  are  A t  <  0.23  s.  To  obtain  statistically  sta¬ 
tionary  solutions  for  these  shear-dominated  flows  re¬ 
quires  about  200  000  time  steps.  Statistics  are  obtained 
by  combined  spatial  and  temporal  averaging  over  the 
last  3-4  h  of  the  simulations.  To  reduce  the  computa- 


Fig.  4.  An  x-z  slice  of  the  conformal  mesh  in  the  lowest  50  m  used  in  the  LES  of  turbulent 
flow  over  water  waves.  The  amplitude  of  the  wave  a  =  1.6  m  and  the  entire  computational 
domain  is  ( xL ,  yL,  Zl)  =  (1200,  1200,  800)  m.  The  cell  aspect  ratio  is  distorted  by  the  plotting 
scales.  In  the  computational  coordinates  ( £,  tj,  f),  surfaces  of  constant  f,  i.e„  £-r)  planes,  follow 
the  underlying  wavy  shape  at  f  =  0  and  smoothly  blend  into  x-y  planes  as  f  increases  away 
from  the  boundary. 
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tional  costs  many  of  the  simulations  are  created  from  an 
initial  seed  run  containing  fully  developed  turbulence. 
The  computational  burden  is  further  increased  by  the 
iterative  pressure  solution  (see  the  appendix),  which 
increases  the  cost  of  the  present  simulations  by  a  factor 
of  2  compared  to  LES  over  a  flat  surface  with  no  re¬ 
solved  surface  undulations. 


4.  LES  results 

Coupling  of  the  wind  fields  to  the  underlying  wavy 
surface  is  found  in  instantaneous  fields  (see  section  4a) 
and  also  in  low-order  statistical  moments  (see  section 
4b).  To  interpret  these  results  it  is  helpful  to  first 
present  the  ensemble  average  equations  of  motion 
above  an  imposed  wave.  The  LES  equations  for  the 
resolved  Cartesian  velocity  components  u  =  (17,  v,  tv), 
where  the  overbar  indicates  spatial  filtering,  are  formu¬ 
lated  in  surface-fitted  wave -following  coordinates  (£,  tj, 
£)  [see  (A12)  in  the  appendix  and  Fig.  4].  In  these  co¬ 
ordinates  and  in  a  frame  of  reference  moving  with  the 
wave  speed  c,  the  ensemble  average  momentum  equa¬ 
tions  for  mean  (77,  v)  in  the  horizontally  homogeneous 
limit  are 


d  /u\  d  /  l, 

-l-r)+-AWu  +  T11^ 


dt  \J 


dC 


Tl3 7+p*J 


=/' 


v~V„ 


(la) 


d  /v\  d 

At\j/  + 


W-  4.  &  4  & 

Wv+  t12j  +  T23y 


(lb) 


where  we  have  assumed  the  special  situation  of  a  2D 
surface  wave.  In  (1),  W  is  the  contravariant  flux  velocity 
normal  to  a  £  surface,  p*  =  plp0  +  2e/3  is  the  pressure, 
r fj  are  subgrid-scale  momentum  fluxes,  e  =  t,,/2  is  the 
subgrid-scale  energy,  pD  is  a  reference  air  density,  /  is 
the  Coriolis  parameter,  (l,x-  £-)  are  grid  metrics,  and  /  is 
the  Jacobian  of  the  grid  transformation.  Also,  angled 
brackets  (•)  denote  a  spatial  average  over  (|,  tj)  coor¬ 
dinates  along  lines  of  constant  £.  For  steady  flow  the 
mean  Ekman  motions  following  a  wavy  surface  result 
from  the  force  balance  between  vertical  divergence  of 
total  (vertical)  momentum  flux  and  geostrophic  pres¬ 
sure  gradients  (— /  Vg,  f  Ug,  0).  The  momentum  flux 
terms  that  appear  on  the  left-hand  side  of  these  equa¬ 
tions  contract  to  their  common  form  for  a  PBL  above  a 
flat  uniform  surface  (e.g.,  Garratt  1992)  when  the  com¬ 
putational  |  gridlines  tend  to  flat  surfaces  in  physical 
space.  As  shown  in  Fig.  4,  near  the  surface  the  mesh 
oscillates  with  the  underlying  wave  but  the  computa¬ 


tional  coordinates  (£,  £)  are  generally  aligned  with  the 
Cartesian  coordinates  ( x ,  z).  Above  z  ~  50  m  the  hori¬ 
zontal  gridlines  are  effectively  level  and  then  — >  xh 
lxIJ  — >  0  so  that  W  -4  vv.  Note  at  £  =  0,  the  wavy  PBL 
Ekman  Eqs.  (1)  contain  an  explicit  contribution  from 
the  pressure  distribution  along  the  wave  p*  'QLJ  —  —p* 
Z?.  This  term  accounts  for  the  resolved  form  stress  (i.e., 
the  drag  or  thrust)  of  the  underlying  wave  (e.g.,  Sulli¬ 
van  et  al.  2000).  As  discussed  below,  the  sign  of  the 
surface  form  stress  depends  on  the  wind-wave  orienta¬ 
tion  and  wave  age. 

We  note  ensemble  statistics  can  also  be  obtained  at 
constant  z  by  first  interpolating  the  wave-following 
computational  results  to  level  z  planes  and  then  aver¬ 
aging.  This  mimics  the  observational  approach  to  gath¬ 
ering  statistics.  Similar  results  are  obtained  from  the 
two  approaches  at  the  same  nominal  z  height  (Sullivan 
et  al.  2000).  However,  results  in  a  constant  z  coordinate 
system  do  not  provide  information  about  flow  dynamics 
between  and  below  wave  crests. 


a.  Flow  visualization 

Extensive  visualization  of  the  LES  solutions  high¬ 
lights  the  impact  of  moving  and  stationary  waves  on 
surface-layer  flow  dynamics,  and  more  broadly,  on  the 
overall  PBL.  The  flow  response  to  stationary  bumps, 
wind  following  waves,  and  wind  opposing  waves  is  radi¬ 
cally  different  as  illustrated  in  Figs.  5  and  6.  Inspection 
of  the  snapshots  (Fig.  5)  shows  an  unexpected  coupling 
between  the  horizontal  winds  and  waves  in  the  situation 
with  wind  following  waves  (wave  age  clUg  «=  2.5).  In 
this  case  a  coherent  pattern  of  accelerated  winds 
greater  than  Ug  occurs  in  the  region  above  each  wave 
trough,  for  5  <  z  <  25  m;  the  fastest  local  winds  occur 
at  z  =  15  m  above  the  mean  water  level  and  the  u  winds 
are  slowest  over  the  wave  crests.  The  formation  and 
coherence  of  a  near-surface  wind  maximum  (or  super- 
geostrophic  jet)  in  this  neutrally  stratified  flow  results 
from  the  coupling  with  fast-moving  surface  waves  and 
not  from  the  interaction  between  turbulence  and  stable 
stratification  as  in  a  nocturnal  boundary  layer  (e.g., 
Saiki  et  al.  2000;  Beare  et  al.  2006).  The  spatial  patterns 
of  the  surface-layer  winds  in  the  presence  of  bumps  or 
waves  opposing  the  wind  are  noticeably  different. 
Overall  the  surface-layer  winds  in  these  two  cases  are 
slower  by  roughly  a  factor  of  4  compared  to  the  situa¬ 
tion  of  wind  following  waves  and  are  always  weaker 
than  the  geostrophic  wind.  The  slowest  surface-layer 
winds  occur  in  the  wave  troughs  in  case  ON  or  on  the 
windward  face  of  the  wave  for  case  BN.  The  organiza¬ 
tion  of  the  resolved  vertical  velocity  w,  shown  in  Fig.  6, 
depends  on  the  wind-wave  orientation  and  wave  age 
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Fig.  5.  Contours  of  the  u  component  of  the  horizontal  wind  field  for  cases  with  moving  and  stationary  surface  waves.  The  nondi- 
mensional  field  shown  is  ulUg.  (top)  Wind  following  waves;  (middle)  wind  opposing  waves;  and  (bottom)  stationary  bumps.  For  each 
case  the  geostrophic  wind  ( Ug ,  Vg)  =  (5,  0)  m  s_1  and  the  wave  slope  ak  =  0.1  where  the  wave  amplitude  a  =  1.6  m.  In  the  top  and 
middle  panels  the  wave  phase  speed  c  =  12.5  m  s-1.  The  color  bar  changes  between  the  top  and  middle  panels.  Note  the  supergeo- 
strophic  winds  near  the  surface  in  the  top  panel. 


similar  to  the  horizontal  wind.  For  wind  following 
waves,  negative  (positive)  patches  of  w  form  upstream 
(downstream)  of  the  wave  crest.  This  pattern  switches 
for  wind  opposing  waves  and  stationary  bumps.  A  com¬ 
parison  of  cases  FN,  ON,  and  BN  shows  wave  propa¬ 
gation  enhances  the  coherence  of  the  vertical  velocity 
and  alters  the  phase  relationship  between  (n.  w)  com¬ 
pared  to  stationary  bumps.  This  implies  fast-moving 


waves  can  impact  the  distribution  of  vertical  momen¬ 
tum  flux  as  discussed  in  section  4b.  We  mention  that  the 
flow  patterns  (not  shown)  in  the  presence  of  high  winds 
Ug  =  12.5  m  s“\  which  are  representative  of  smaller 
wave  age,  are  qualitatively  similar  to  results  for  flow 
over  stationary  bumps  (Sullivan  et  al.  2000). 

Our  interpretation  of  the  cause  of  the  u-w  flow  pat¬ 
terns  in  the  surface  layer  is  based  on  the  structure  of  the 
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Fig.  6.  Contours  of  the  w  component  of  the  vertical  wind  field  for  cases  with  moving  and  stationary  surface  waves  as  in  Fig.  5.  The 
nondimensional  field  shown  is  wlUg.  (top)  Wind  following  waves;  (middle)  wind  opposing  waves;  and  (bottom)  stationary  bumps. 


near-surface  pressure  field  p*  shown  in  Figs.  7  and  8, 
and  in  particular  the  sign  of  the  form  stress  (i.e.,  the 
surface  drag)  induced  by  the  waves.  In  these  plots  we 
show  the  phase-averaged  pressure  signal  [p*(x,  z)]  =  Jv 
p*{x,  y,  z)dylyL  normalized  by  Uj,r  In  the  simulations 


2  Our  choice  of  normalization  based  on  Ug  instead  of  friction 
velocity  h*  results  from  the  observation  that  depending  on  wind- 
wave  alignment  and  wave  age  the  vertical  momentum  flux  can 
change  sign  or  vary  appreciably  with  z  in  the  surface  layer,  which 
leads  to  a  poorly  defined  u *. 


a  coherent  pattern  of  positive  and  negative  pressure 
correlated  with  the  wave  crests  and  troughs  develops 
and  extends  well  above  the  surface.  Comparison  of  the 
three  cases  shows  1)  the  weakest  pressure  fluctuations 
occur  in  the  case  with  wind  following  the  waves;  2)  wind 
opposing  waves  generates  the  most  vigorous  fluctua¬ 
tions,  which  can  extend  to  a  height  of  27rz/A  >  2;  and 
3)  there  is  a  subtle  asymmetry  in  the  pressure  field 
relative  to  the  underlying  wave  depending  on  the  wave 
age  and  wind-wave  orientation  that  leads  to  the  form 
stress.  In  case  FN,  the  negative  pressure  pattern  is 
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Fig.  7.  Contours  of  the  nondimensional  and  y-averaged  pressure  field  {p*)IU2g  close  to  the 
water  surface  for  cases  with  moving  and  stationary  waves.  The  winds  are  from  left  to  right. 
Negative  contours  are  indicated  by  dashed  lines,  (top)  Wind  following  waves;  (middle)  wind 
opposing  waves;  and  (bottom)  stationary  bumps.  The  vertical  and  horizontal  coordinates  are 
made  dimensionless  with  the  surface  wavelength  A. 


shifted  slightly  behind  the  wave  crest  (x/ A  <  1);  hence, 
the  integration  of  the  surface  pressure  over  the  wave 
acts  in  the  positive  x  direction,  implying  a  thrust  on  the 
winds.  Meanwhile  in  cases  ON  and  BN,  the  negative 
pressure  minimum  is  shifted  slightly  ahead  of  the  wave 
crest  (x/A  >  1)  and  then  the  surface  form  stress  acts  as 
a  drag  on  the  surface  winds  as  expected  for  stationary 
roughness.  The  magnitude  and  sign  of  the  form  stress 
reflects  the  change  in  character  of  the  surface-layer  tur¬ 
bulence. 

The  pattern  of  surface  pressure  in  the  case  with  wind 
following  waves  observed  here  in  LES  of  a  full  PBL  is 
qualitatively  similar  to  the  predictions  from  linear 
analysis  (Belcher  and  Hunt  1998;  Cohen  and  Belcher 
1999),  second-order  closure  (Gent  1977;  Kudryavtsev 
and  Makin  2004),  and  direct  numerical  simulations 
(Sullivan  et  al.  2000).  All  predict  that  for  large  values  of 
wave  age  the  form  stress  acts  as  a  thrust  on  the  winds; 


Fig.  8.  Streamwise  x  variation  of  the  nondimensional  and  y- 
averaged  surface  pressure  for  cases  with  moving  and  stationary 
surface  waves,  (bottom)  The  underlying  wave.  Here,  diamonds 
indicate  wind  following  waves,  triangles  indicate  wind  opposing 
waves,  and  squares  indicate  stationary  bumps.  The  horizontal  co¬ 
ordinate  is  made  dimensionless  by  the  wavelength  A. 
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Fig.  9.  Snapshot  of  the  resolved  vertical  momentum  flux  u'w'IU2  in  an  x-y  plane  at  z  «=  15  m  above  the  surface,  (a)  Flat  z0  surface, 
(b)  stationary  bumps,  (c)  wind  following  waves,  and  (d)  wind  opposing  waves.  The  wave  and  wind  conditions  are  as  described  in 
Fig.  5. 


this  results  from  an  asymmetrical  pressure  distribution 
with  the  minimum  negative  pressure  forward  of  the 
wave  crest. 

Surface  waves  impact  the  instantaneous  velocity 
and  pressure  fields,  and  thus  we  next  examine  how 
surface  waves  modulate  the  important  momentum- 
flux-carrying  coherent  structures  in  the  surface  layer. 
The  flow  visualization  in  Fig.  9  compares  the  instanta¬ 
neous  (resolved)  momentum  flux  u'w'  at  a  nominal 
height  of  z  =  15  m  (or  £/z,-  =  0.0375)  above  different 


surfaces.3  At  this  z  the  horizontal  £;  gridlines  are  effec¬ 
tively  level,  and  the  momentum  flux  is  dominated  by 
resolved  fluctuations.  We  observe  over  a  flat  za  bound¬ 
ary  the  bulk  of  the  negatively  signed  vertical  momen¬ 
tum  flux  is  carried  by  a  few  sparsely  distributed  struc¬ 
tures  aligned  in  the  mean  wind  direction.  Similar  elon- 


3  Here  ()'  denotes  a  deviation  from  a  horizontal  average,  that  is, 
a  turbulent  fluctuation. 
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gated  flux-carrying  structures  are  observed  in  direct 
numerical  simulations  over  a  smooth  wall  (e.g.,  Adrian 
and  Liu  2002),  in  other  LES  (e.g.,  Lin  et  al.  1996; 
Moeng  and  Sullivan  1994),  and  also  in  outdoor  obser¬ 
vations  (e.g.,  Hommema  and  Adrian  2003).  Fast- 
moving  waves  leading  or  opposing  the  wind  destroy  the 
coherence  of  these  streaky  near-wall  structures.  For 
wind  following  waves,  the  momentum  flux  structures  in 
the  surface  layer  are  weak  and  carry  slightly  positive 
flux  and  impact  the  ensemble  average  profile  (shown 
later  in  Fig.  11).  The  scale  of  the  structures  is  observed 
to  be  linked  to  the  horizontal  scale  of  the  waves.  Chang¬ 
ing  the  direction  of  wave  propagation  relative  to  the 
winds  drastically  alters  the  momentum  flux  patterns. 
Turbulent  structures  carrying  large  amounts  of  positive 
and  negative  momentum  strongly  correlated  with  the 
motion  of  the  underlying  waves  are  observed  in  Fig.  9d. 
Additional  visualization  shows  that  Ti'w'  induced  by  op¬ 
posing  waves  remains  coherent  well  above  the  surface 
layer  and  appears  to  interact  with  the  background  PBL 
turbulence.  The  structural  features  of  the  momentum 
flux  in  case  ON  are  consistent  with  the  velocity  and 
pressure  fields  discussed  previously.  Wave  age  and 
wind-wave  orientation  are  then  clearly  important  for 
momentum  flux  generation  since  stationary  bumps  of 
the  same  amplitude  as  the  moving  waves  considered 
here  generate  flux  structures  more  comparable  to  those 
above  a  flat  surface. 

Animations  of  these  and  other  LES  solutions  dem¬ 
onstrate  that  the  structure  of  the  velocity,  pressure,  and 
momentum  flux  fields  are  persistent  in  time  and  robust 
to  reductions  in  surface  roughness  za  and  the  presence 
of  slight  surface  heating.  As  we  illustrate  later  the  im¬ 
pact  of  surface  waves,  especially  in  the  case  of  wind 
following  waves,  is  not  confined  to  the  surface  layer  but 
can  extend  over  the  PBL.  In  this  case  the  mean  shear  is 
weak  between  the  top  of  the  surface-layer  jet  and  zt, 
which  reduces  turbulence  production  in  the  bulk  of  the 
PBL.  Meanwhile  the  same  wave  moving  in  opposition 
to  the  wind  acts  as  a  large  drag  element  slowing  the 
surface-layer  winds  and  generating  vigorous  turbulence 
that  fills  the  PBL.  The  flow  patterns  found  here  in  the 
presence  of  moving  waves  are  in  contrast  to  flow  over  a 
stationary  hill  (Belcher  and  Hunt  1998)  and  suggest 
propagating  waves  can  in  certain  circumstances  modify 
the  overlying  turbulent  flow  over  the  bulk  of  what  is 
traditionally  referred  to  as  the  PBL  surface  layer,  cor¬ 
responding  to  approximately  zb<  z  <  0.1  z,-. 

b.  Vertical  profiles  of  winds,  momentum  fluxes,  and 
variances 

Vertical  profiles  of  the  time-  and  space-averaged 
mean  winds  above  a  flat  za  surface,  stationary  bumps, 


Fig.  10.  Vertical  profiles  of  the  horizontal  components  of  the 
mean  wind  («),  (v)  for  flow  over  waves.  The  spatial  averaging  is 
carried  out  along  constant  f  surfaces,  where  f  is  the  mean  height 
above  the  wave.  The  nominal  boundary  layer  depth  z;  =  400  m. 
The  cases  are  as  follows:  dotted  line,  no  waves;  squares,  stationary 
bumps;  diamonds,  wind  following  waves;  triangles,  wind  opposing 
waves;  and  dashed  line,  slight  convection  with  wind  following 
waves. 


and  moving  waves  are  compared  in  Fig.  10.  As  antici¬ 
pated  based  on  the  flow  visualization  in  section  4a,  the 
(u,  d)  wind  profiles  above  stationary  bumps  and  wind 
opposing  waves  are  broadly  similar  to  those  above  a  flat 
zQ  surface;  the  mean  u  profile  is  positively  sheared  over 
the  entire  PBL  and  stationary  bumps  and  waves  oppos¬ 
ing  the  wind  generate  large  mean  vertical  gradients.  In 
the  situation  of  wind  following  waves  the  structure  of 
the  mean  wind  profiles  is  radically  different;  the  u  wind 
profile  exhibits  a  low-level  jet  (u)/Ug  >  1.1  near  z  ~  20 
m  and  the  sign  of  the  v  profile  is  switched  compared  to 
a  flat  surface.  Above  the  low-level  maximum  the  winds 
smoothly  transition  to  the  geostrophic  wind  with  (»)  «* 
Ug  for  z  &  200  m.  With  small  surface  heating  the  low- 
level  jet  nearly  disappears,  (u)  Ug ,  and  d(u)ldz  553  0 
over  the  bulk  of  the  PBL,  for  10  m  <  z  <  400  m.  Hence 
in  both  the  FN  and  FC  cases  the  mean  shear  above  the 
jet  is  either  slightly  negative  or  nearly  zero.  Below  the 
height  of  the  low-level  maximum  the  u  winds  decrease 
sharply  in  order  to  match  the  surface  boundary  condi¬ 
tions.  Because  of  the  interaction  with  the  wave  field  the 
horizontal  wind  direction  depends  on  wave  state.  Given 
the  orientation  of  the  geostrophic  wind,  Ug  parallel  to 
the  x  direction,  the  surface  winds  at  the  standard  ref¬ 
erence  height  z  ~  10  m  turn  to  the  left  as  expected  for 
flow  over  a  flat  za  surface  as  (v)  >  0.  The  degree  of 
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Fig.  11.  Vertical  profiles  of  the  nondimensional  vertical  mo¬ 
mentum  flux  for  flow  over  waves,  (a)  The  sum  (Wu  +  r, , £x U  + 
t 13f,//  +  p*£xIJ)  X  100 lUg,  and  (b)  the  sum  (Wv  +  r12^//  + 
t 23izU)  x  100 /Ug.  The  spatial  averaging  is  carried  out  along  con¬ 
stant  f  surfaces,  where  f  is  the  mean  height  above  the  wave  and  z,. 
=  400  m.  The  cases  are  as  follows:  dotted  line,  no  waves;  squares, 
stationary  bumps;  diamonds,  wind  following  waves;  triangles, 
wind  opposing  waves;  and  dashed  line,  slight  convection  with  wind 
following  waves. 


turning  increases  for  stationary  bumps  consistent  with 
their  larger  surface  form  stress.  However,  an  opposite 
trend  is  observed  in  the  presence  of  wind  following 
waves;  the  winds  turn  slightly  to  the  right  in  the  surface 
layer  and  (v)  <  0,  so  that  they  are  nearly  aligned  with 
the  wave  propagation  direction.  The  rightward  turning 
of  the  wind  is  a  consequence  of  the  Ekman  balance  for 
momentum  above  waves  discussed  below.  This  LES 
prediction  is  at  least  qualitatively  similar  to  the  obser¬ 
vations  in  light  winds  above  swell  reported  by  Grachev 
et  al.  (2003). 

Surface  waves  modify  the  momentum  balance  in  the 
atmospheric  PBL  as  shown  in  Figs.  11  and  12.  In  Fig.  11 
the  vertical  distribution  of  the  sum  of  flux  contributions 
on  the  left  hand  side  of  the  Ekman  Eq.  (1)  is  shown:  in 
order,  these  terms  are  resolved  momentum  flux,  sub- 
grid-scale  contributions,  and  pressure  stress.  In  the 
cases  with  stationary  bumps  and  wind  opposing  waves 
the  distribution  of  the  two  components  of  the  vertical 
momentum  flux  are  as  expected  for  a  turbulent  PBL 
above  a  rough  surface;  the  dominant  u  momentum  flux 
is  negative  with  positive  vertical  divergence.  Fast- 
moving  waves  leading  the  wind  greatly  alter  the  mo¬ 
mentum  flux  distribution;  the  u  component  is  slightly 
positive  with  negative  vertical  divergence  while  the  sign 


Fig.  12.  Vertical  profiles  of  the  normalized  pressure  stress 
(p*£xIJ)  X  100/f/g  near  the  water  surface  for  flow  over  waves. 
The  pressure  is  averaged  along  constant  f  surfaces  with  zt  =  400 
m.  The  cases  are  as  follows:  squares,  stationary  bumps;  diamonds, 
wind  following  waves;  triangles,  wind  opposing  waves;  and  dashed 
line,  slight  convection  with  wind  following  waves.  Note  for  cases 
with  wind  following  waves  the  waves  impart  a  forward  (positive) 
thrust  on  the  winds. 


and  vertical  divergence  of  the  v  component  are  oppo¬ 
site  to  their  counterparts  in  a  flat  PBL.  The  cause  of  this 
unexpected  behavior  can  be  traced  to  the  pressure 
stress  variation  shown  in  Fig.  12.  At  the  wave  surface, 
fast-moving  waves  impart  a  positive  forward  thrust  on 
the  winds  opposite  to  that  in  flow  over  stationary 
bumps.  In  other  words  there  is  significant  momentum 
transfer  from  the  ocean  to  the  atmosphere.  The  surface 
thrust  from  the  waves  is  a  large  component  of  the  mo¬ 
mentum  flux  balance  and  acts  counter  to  the  usual  drag 
induced  by  surface-generated  turbulence.  The  vertical 
distribution  of  pressure  stress  above  the  surface  f  >  0  is 
a  consequence  of  formulating  the  Ekman  flux  budget 
(1)  in  wave-following  coordinates.  Its  smooth  mono¬ 
tonic  variation  with  height  shows  that  the  surface  asym¬ 
metry  of  the  pressure  contours  with  respect  to  the  un¬ 
derlying  wave  field  (see  Fig.  7)  persists  with  increasing 
distance  z  from  the  surface. 

The  variation  and  signs  of  both  components  of  the 
momentum  flux  and  mean  winds  are  consistent  with  the 
formation  of  a  low-level  jet  and  are  mandatory  in  order 
to  achieve  a  steady  balance  between  the  pressure  gra¬ 
dient  forcing  and  momentum  flux  divergence.  With 
wind  following  waves  the  Ekman  balance  of  terms  is 
opposite  to  that  of  a  conventional  PBL,  the  stress  di¬ 
vergence  serves  to  accelerate  the  u  component  of  the 
wind  while  the  pressure  gradient  acts  to  retard  the  flow. 
Finally,  notice  that  with  small  amounts  of  convection 
the  vertical  //-momentum  flux  is  small  but  clearly  posi¬ 
tive  over  the  vertical  extent  30  m  <  z  <  z,- .  This  is  in 
contrast  to  a  PBL  over  a  land  surface  driven  by  shear 
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Fig.  13.  Vertical  profiles  of  the  nondimensional 
resolved  variance  components  (uj  u'i)IU2g  (no  sum 
over  /).  The  results  are  shown  in  linear-logarith¬ 
mic  coordinates  with  the  vertical  axis  nondimen- 
sionalized  by  the  initial  height  of  the  inversion  zt 
=  400  m.  The  cases  are  as  follows:  dotted  line,  no 
waves;  squares,  stationary  bumps;  diamonds,  wind 
following  waves;  triangles,  wind  opposing  waves; 
and  dashed  line,  slight  convection  with  wind  fol¬ 
lowing  waves.  The  horizontal  arrow  shows  the 
vertical  location  of  the  low-level  wind  maximum 
in  case  FN. 


and  convection  where  the  momentum  flux  in  the  upper 
PBL  is  negative  (e.g.,  Moeng  and  Sullivan  1994).  We 
speculate  surface  convection  transports  positive  signed 
vertical  momentum,  generated  by  the  wave  field,  to  the 
upper  regions  of  the  PBL. 

For  a  horizontally  homogeneous  PBL  above  a  flat 
surface,  the  turbulent  kinetic  energy  (TKE)  budget 
(e.g.,  Moeng  and  Wyngaard  1989;  Moeng  and  Sullivan 
1994)  contains  two  main  sources  of  energy,  namely, 
shear  production  T  =  — (u'w')  •  d(u)ldz  and  buoyancy 
S  =  g/0o  (w'Q').4  Here  (P>  0  since  (u'w')  and  d(u)/dz  are 
generally  opposite  in  sign.  In  our  shear-dominated 


4  In  the  definitions  of  S’  and  $  the  ()'  denotes  a  turbulent  fluc¬ 
tuation  and  the  fluctuating  velocity  vector  ii'  =  («',  v',  w'). 


PBLs  the  presence  of  surface  waves  alters  the  shear 
production  mechanism  and  hence  TKE.  These  changes 
are  reflected  in  the  (resolved)  component  variances 
( u '  u',  v'  v ',  w'  w')  shown  in  Fig.  13.  Wave  influences 
dominate  near  the  surface  but  also  impact  the  distribu¬ 
tion  of  turbulence  energy  over  the  bulk  of  the  PBL.  In 
the  neutral  case,  with  wind  following  waves  a  near¬ 
surface  velocity  maximum  is  generated  with  slightly  su- 
pergeostrophic  winds  (u)IUg  1.1  (see  Fig.  10).  As  a 
result  the  shear  between  the  surface  wind  maximum 
and  the  PBL  top  is  near  zero.  Coupled  with  small  ver¬ 
tical  momentum  fluxes  (see  Fig.  11),  the  shear  produc¬ 
tion  ?  is  minimal  over  the  bulk  of  the  PBL.  Note  in  Fig. 
13  when  £/z,-  >0.1  the  smallest  variances  occur  in  case 
FN.  Swell  propagating  in  the  wind  direction  then  has  a 
significant  impact  on  the  turbulence  level  in  the  neutral 
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PBL  as  the  modification  of  the  turbulence  production 
mechanism  in  the  surface  layer  leads  to  a  turbulence 
collapse  in  the  overall  PBL.  Surface  convection,  how¬ 
ever,  still  generates  significant  TKE  in  the  PBL  in  the 
presence  of  waves  as  shown  in  Fig.  13.  Stationary 
bumps  and  waves  opposing  the  surface  wind  both  gen¬ 
erate  turbulence  variances  larger  than  a  neutrally  strati¬ 
fied  flat  z0  surface  consistent  with  their  larger  surface 
form  drag  and  sheared  mean  wind  profiles.  Near  the 
surface,  £/z;-  <  0.1,  and  the  («,  vv )  variances  in  the  pres¬ 
ence  of  moving  waves  are  large  due  to  the  significant 
(irrotational)  motion  of  the  underlying  wave  field. 

These  LES  predictions  in  the  marine  surface  layer 
are  qualitatively  supported  by  the  observations  of 
Smedman  et  al.  (1999)  who  find  that  turbulence  pro¬ 
duction  is  significantly  reduced  in  the  presence  of  wind 
following  waves.  Thus  TKE  in  the  marine  PBL  depends 
on  wind-wave  orientation,  wave  age,  and  generally  on 
the  structure  of  the  wave  field. 

5.  Momentum  fluxes  from  CBLAST  and  LES 

Our  LES  results  predict  that  the  winds,  turbulence 
fluxes,  and  variances  as  well  as  their  mean  profiles  de¬ 
pend  on  bulk  properties  of  the  wave  field,  that  is,  wave 
age  and  wind-wave  orientation.  These  computational 
results  provide  motivation  to  search  for  wave  influences 
in  measured  wind  fields  from  the  CBLAST  field  cam¬ 
paign.  Compared  to  real  seas,  the  wave  fields  in  the 
LES  are  highly  idealized;  for  example,  they  do  not  in¬ 
clude  multicomponents,  three-dimensionality,  and 
time-varying  wave  amplitudes  and  phases.  Hence,  we 
expect  wave  influences  to  be  more  subtle  and  difficult 
to  isolate  in  observations. 

First  we  interrogate  the  CBLAST  database  searching 
for  cases  with  winds  and  waves  that  conform  to  the  LES 
idealizations  for  more  detailed  analysis.  Based  on  a  cri¬ 
terion  of  wind-wave  angle  —30°  <  <f>  <  30°,  approxi¬ 
mately  100  periods  of  60  min  in  duration  are  identified 
as  cases  of  wind  following  waves.  Unfortunately,  nu¬ 
merous  clean  cases  with  waves  directly  opposing  the 
winds  are  not  present  because  of  possible  flow  distor¬ 
tion  from  the  ASIT  superstructure  (see  Fig.  1).  By  ex¬ 
panding  the  wind-wave  angle  to  130°  <  <j>  <  230°  a 
limited  number  of  cases  are  identified  as  wind  opposing 
waves — just  18  periods  of  20-min  duration.  In  the  data 
screening,  no  limits  are  placed  on  the  range  of  atmo¬ 
spheric  stratification,  but  we  note  that  stable  stratifica¬ 
tion  is  a  potential  source  of  variability  (Smedman  et  al. 
1997).  For  the  selected  cases  with  wind  following  waves 
shown  in  Fig.  14,  the  wave  age  CpIUacos  <fi  spans  a  large 
range,  approximately  [1,  8].  This  subset  of  the  CBLAST 
data  is  dominated  by  fast-moving  swell  (or  old  seas) 


Fig.  14.  Frequency  of  wave  age  for  selected  cases  with  wind 
following  waves.  The  solid  vertical  bars  show  the  frequency  of 
occurrence  and  the  solid  line  is  the  cumulative  probability  sum 
1  —  fo  p(x')dx’ ,  where  p(x)  is  the  probability  density  function. 

with  approximately  a  75%  probability  of  wave  age 
greater  than  wind-wave  equilibrium;  the  high  probabil¬ 
ity  of  old  seas  in  this  subset  of  data  is  similar  to  that  for 
the  entire  CBLAST  dataset.  Hence  we  expect  wave 
influences  are  present  in  this  subset  of  the  CBLAST 
data. 

LES  predicts  the  surface-atmosphere  momentum  ex¬ 
change  depends  on  wave  state.  To  expose  and  quantify 
this  dependence  in  the  CBLAST  data  a  quadrant  analy¬ 
sis  of  the  observed  vertical  momentum  flux  is  per¬ 
formed.  This  conditional  sampling  technique,  first  used 
with  observational  data  by  Chambers  and  Antonia 
(1981)  and  later  by  Smedman  et  al.  (1999),  separates 
the  turbulent  momentum  flux  u'w'  into  four  categories 
(quadrants)  according  to  the  sign  of  the  two  fluctuating 
velocity  components  as  sketched  in  Fig.  15.  In  the  sur¬ 
face  layer  of  a  rough  wall  boundary  layer  the  net  (av¬ 
erage)  momentum  flux  {u'w')  <  0  and  is  dominated  by 
sweeps  and  ejections  associated  with  motions  in  quad¬ 
rants  Q2  and  Q4.  Positive  flux  contributions  from  the 
interaction  quadrants  Q1  and  Q3  are  less  frequent  and 
weaker  in  magnitude. 

A  quadrant  analysis  of  the  vertical  momentum  flux 
obtained  in  CBLAST  and  from  our  idealized  LES  data 
is  displayed  in  Fig.  16.  The  CBLAST  results  for  wind 
following  waves  are  averages  over  the  four  vertical 
sonic  positions  z  =  [4.0,  6.5, 10.0,  18.0]  m  and  wave  age 
bins  of  width  equal  to  0.09.  In  the  cases  with  wind  op¬ 
posing  waves,  the  observational  results  are  only  aver¬ 
aged  over  the  four  sonic  positions  owing  to  the  limited 
dataset.  For  comparison  we  also  display  observational 
results  for  flow  over  stationary  (terrestrial)  roughness 
(Sullivan  et  al.  2003)  and  from  Smedman  et  al.  (1999) 
obtained  at  an  independent  marine  site.  The  results  are 
presented  in  terms  of  the  normalized  ratio  of  negative 
to  positive  momentum  flux  quadrants  Qr  =  —  (Q2  + 
Q4)/(Q1  +  Q3)  for  varying  wave  age  CpIUa  cos  4>.  We 
find  the  quadrant  ratio  Qr  to  be  a  robust  statistical  mea- 
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Fig.  15.  Decomposition  of  the  vertical  momentum  flux  into 
quadrants  (Ql,  Q2,  Q3,  Q4)  based  on  the  sign  of  the  fluctuating 
horizontal  and  vertical  velocity  ( u ' ,  w'). 


sure  that  exposes  the  nature  of  the  underlying  surface, 
and  in  the  present  analysis  brings  out  the  wave  depen¬ 
dence.  The  CBLAST  results  contain  scatter  but  the 
quadrant  flux  ratio  clearly  contains  wave  influences,  a 
distinct  downward  trend  for  increasing  wave  age  CpIUa 
cos  <f>  >  1.  Our  interpretation,  based  on  the  LES  results, 
is  that  under  low  winds  the  fast-moving  components  of 
the  wave  field  enhance  the  upward  (positive)  momen¬ 
tum  transport  from  the  ocean  to  the  atmosphere  and 
this  momentum  appears  in  the  positively  signed  flux 
quadrants  (Ql,  Q3).  At  a  sufficiently  large  wave  age  a 
near  balance  between  negative  and  positive  flux  contri¬ 
butions  is  reached,  implying  zero  surface  drag.  The 
quadrant  momentum  flux  distributions  are  a  conse¬ 
quence  of  competing  effects;  fast-moving  waves  gener¬ 
ate  positive  momentum  flux  while  small  slow-moving 
waves  act  similar  to  conventional  roughness  elements. 
Also  the  effects  of  fast-moving  waves  on  momentum 
transport  are  not  confined  to  the  first  measurement 
level,  z  =  4.0  m,  but  extend  over  the  bulk  of  the  surface 
layer,  up  to  at  least  z  =  18.0  m,  in  agreement  with  the 
LES.  The  few  observations  reported  by  Smedman  et  al. 
(1999)  also  follow  a  similar  trend  with  wave  age  as  the 
CBLAST  results.  Notice  Qr  appears  to  asymptotically 
approach  a  value  measured  at  a  rough  land  site  for 
wave  age  approaching  zero,  that  is,  a  field  of  young 
developing  waves  generates  a  distribution  of  momen¬ 
tum  flux  broadly  similar  to  stationary  roughness.  Al¬ 
though  Qr  provides  information  as  to  the  distribution  of 
momentum  flux  it  does  not  provide  scale  information. 


Fig.  16.  Quadrant  analysis  of  the  vertical  momentum  flux  in  the 
marine  surface  layer  for  varying  wave  age  with  wind  following  and 
opposing  waves.  CBLAST  results  are  indicated  by  crossed  circles. 
For  comparison  we  show  observations  of  Smedman  et  al.  (1999), 
denoted  by  X,  and  results  for  flow  over  stationary  roughness  (note 
wave  age  =  0)  (Sullivan  et  al.  2003),  indicated  by  an  open  square 
with  an  error  bar.  LES  results  at  z  »»  15  m  above  the  surface  are 
indicated  by  large  filled  symbols:  circles,  stationary  bumps;  dia¬ 
monds,  wind  opposing  waves  (note  wave  age  <0);  left-pointing 
triangles,  wind  following  waves;  right-pointing  triangles,  wind  fol¬ 
lowing  waves  plus  convection;  and  triangles,  wind  following  waves 
for  very  light  winds  with  Ug  =  2  m  s  A 

Spectral  analysis  of  surface-layer  winds  in  the  presence 
of  waves  (Drennan  et  al.  1999;  Smedman  et  al.  2003) 
show  that  the  low  wavenumbers  (or  frequencies)  are 
modified  by  swell,  that  is,  for  wind  following  waves. 

The  LES  predictions  for  the  distribution  of  vertical 
momentum  flux  are  in  general  good  agreement  with  the 
observational  trends.  They  mimic  the  observed  varia¬ 
tion  with  wave  age  but  likely  overemphasize  the  wave- 
driven  wind  effects  due  to  the  highly  idealized  and  per¬ 
sistent  nature  of  the  surface  wave  field.  Thus  LES  pre¬ 
dicts  lower  values  of  Qr.  Also,  LES  hints  at  an 
intriguing  wave  effect  on  the  momentum  flux  for  flows 
with  wind  opposing  waves.  Case  ON  with  wave  age 
—4.7  is  characterized  by  high  CD  (see  Tabic  1),  small 
surface  wind  speed,  large  momentum  flux,  and  high 
variance.  Opposing  waves  are  efficient  generators  of 
fluctuation  amplitude,  which  modulates  the  flux¬ 
carrying  structures  compared  to  a  flat  surface  as  shown 
in  Fig.  9.  This  leads  to  a  reduced  value  of  Qr  in  Fig.  16. 
The  LES  predictions  for  wind  opposing  waves  are  also 
verified  by  a  limited  number  of  CBLAST  observations. 
The  present  LES  results  are  further  supported  by  inde¬ 
pendent  direct  numerical  simulations  of  a  wavy  Couette 
flow  (Sullivan  et  al.  2000). 

Finally  the  wave  influences  observed  in  the  vertical 
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Fig.  17.  The  variation  of  the  neutral  drag  coefficient  with  wind  speed  for  wind  following 
waves  in  CBLAST.  The  wave  age  for  these  observations  is  mostly  greater  than  1.2  as  shown 
in  Fig.  14.  The  vertical  locations  are  nominally  z  =  (4.0,  squares;  6.5,  diamonds;  10.0,  circles) 
m.  The  TOGA  COARE  3.0  parameterization  is  the  solid  line. 


momentum  flux  naturally  appear  in  the  bulk  measure¬ 
ment  of  the  sea  surface  drag.  Figure  17  shows  the  varia¬ 
tion  of  the  drag  coefficient  obtained  in  CBLAST  for  all 
cases  with  wind  following  waves.  Again  we  emphasize 
that  these  results  cover  a  range  of  sea  state  but  are 
dominated  by  wave-age  conditions  greater  than  wind- 
wave  equilibrium.  Notice  the  majority  of  the  CD  values 
fall  well  below  the  standard  TOGA  COARE  param¬ 
eterization  especially  at  low  wind  speed.  This  effect  is 
due  to  the  presence  of  the  underlying  swell,  which  in¬ 
duces  upward  momentum  transport  from  the  ocean  to 
the  atmosphere  as  predicted  by  LES.  In  these  cases  the 
wave  field  alters  the  usual  turbulence  production 
mechanism  in  the  marine  surface  layer  and  lowers  the 
drag  coefficient.  Often  the  measured  CD  is  only  50%  of 
the  standard  parameterization  value  and  can  clearly 
approach  negative  values.  The  values  of  CD  from  the 
LES  at  z  =  15  m,  listed  in  Table  1,  are  at  least  quali¬ 
tatively  similar  to  the  CBLAST  measurements.  Com¬ 
pared  to  a  flat  zD-surface,  fast-moving  swell  leading  the 
surface  wind  leads  to  CD  <  0,  while  in  the  presence  of 
opposing  swell,  CD  increases  by  more  than  a  factor  of  4 
in  agreement  with  the  observations  of  Donelan  et  al. 
(1997). 

6.  Conclusions 

Recent  measurements  from  the  Coupled  Boundary 
Layers  Air-Sea  Transfer  (CBLAST)  field  campaign 
(Edson  et  al.  2007)  show  that  the  winds  and  waves  in 
the  marine  surface  layer  are  frequently  in  a  state  of 
disequilibrium  in  light  to  moderate  wind  conditions 
where  Ua  <  10  m  s-1.  Long-wavelength,  fast-moving 
waves  generated  by  distant  storms  often  dominate  the 


local  wave-height  variance  and  spectrum  and  propagate 
in  arbitrary  directions  relative  to  the  local  wind.  In 
terms  of  a  bulk  wave  age  CpIUa  cos  <f>,  where  C„  is  the 
phase  speed  of  the  peak  in  the  wave-height  spectrum 
and  4>  the  wind-wave  angle,  the  wave  age  is  most  often 
either  negative  or  greater  than  the  equilibrium  value  of 
1.2.  In  low-wind  conditions  swell  is  then  an  important 
source  of  variability  in  measurements  of  the  surface 
drag  coefficient  CD. 

To  examine  the  interaction  between  atmospheric  tur¬ 
bulence  and  swell,  a  large-eddy  simulation  (LES) 
model  of  the  planetary  boundary  layer  (PBL)  is  devel¬ 
oped  with  the  capability  of  imposing  propagating  sinu¬ 
soidal  modes  at  its  lower  boundary.  The  code  is  used  to 
simulate  a  variety  of  PBLs  with  an  emphasis  on  situa¬ 
tions  with  wind  following  waves,  wind  opposing  waves, 
and  stationary  bumps.  The  LES  results  illustrate  the 
importance  of  wave  phase  speed  relative  to  wind  speed 
and  the  orientation  of  winds  and  waves.  Surface-layer 
winds  are  modulated  by  the  structure  of  the  near¬ 
surface  pressure  field  (i.e.,  the  resolved  surface  form 
stress).  In  flow  over  stationary  bumps  or  wind  opposing 
waves,  the  resolved  form  stress  is  negative,  while  for 
wind  following  waves,  the  resolved  form  stress  is  posi¬ 
tive.  In  the  latter  situation  LES  predicts  momentum 
transfer  from  the  ocean  to  the  atmosphere  and  the  gen¬ 
eration  of  a  low-level  jet;  the  magnitude  of  the  winds  at 
z  ~  [10,  20]  m  are  about  10%  greater  than  the  geo- 
strophic  wind  and  vary  with  surface  heating.  Our  inter¬ 
pretation  suggests  that  the  jet  formation  results  from  a 
wave-induced  turbulent  momentum  flux  divergence 
that  accelerates  the  flow  and  a  retarding  pressure  gra¬ 
dient,  both  of  which  are  opposite  to  the  momentum 
balance  in  classical  shear  boundary  layers.  In  a  neu- 
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trally  stratified  PBL,  the  presence  of  a  low-level  jet 
reduces  the  mean  shear  between  the  surface  layer  and 
the  PBL  top,  leading  to  a  near  collapse  of  turbulence  in 
the  PBL.  The  mean  wind  profile,  turbulence  variances, 
and  vertical  momentum  flux  are  then  dependent  on  the 
nature  of  the  wave  field,  the  wind-wave  orientation, 
and  wave  age.  The  LES  predictions  for  the  dependence 
of  vertical  momentum  flux  on  wave  age  are  also  found 
in  the  CBLAST  observations.  The  LES  results  with 
moving  waves  show  important  differences  compared 
with  rough-wall  boundary  layers  and  flow  over  station¬ 
ary  bumps  (i.e.,  hills).  Bulk  parameterizations  of  the 
surface  drag  need  to  account  for  wave  state. 
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APPENDIX 

LES  Model 

a.  LES  equations  in  wave-following  coordinates 

An  LES  code  with  the  capability  of  resolving  a  mov¬ 
ing  sinusoidal  mode  imposed  at  its  lower  boundary  was 
developed.  The  computational  approach  is  similar  to 
that  described  in  our  direct  numerical  simulation 
(DNS)  code  (Sullivan  and  McWilliams  2002;  Sullivan  et 
al.  2000):  first  we  translate  the  grid  horizontally  to  take 
out  the  movement  of  the  underlying  surface  and  then 
we  apply  a  grid  transformation  to  the  flow  equations 
mapping  the  physical  domain  into  a  flat  computational 
space  (e.g.,  Anderson  et  al.  1984).  As  is  standard  com¬ 
putational  practice,  the  mapping  is  applied  only  be¬ 
tween  Cartesian  and  computational  coordinates  Xj  — »  f,. 
For  the  atmospheric  PBL  the  working  flow  model  is 
assumed  to  be  unsteady,  3D,  and  described  by  incom¬ 
pressible  Boussinesq  equations  with  large-scale  pres¬ 
sure  gradients  provided  by  geostrophic  winds.  The  gov¬ 
erning  set  of  LES  model  equations  in  Cartesian  coor¬ 
dinates  for  this  flow  is  given  by  Moeng  (1984).  They 
include  transport  equations  for  resolved-scale  (or  spa¬ 
tially  filtered)  velocity  z7,  and  virtual  potential  tempera¬ 
ture  0: 


dllt 

d 

e  ge 

dp-* 

1  dtP 

dt  ~ 

-  fo.  ( UjUt  +  T ij) 

8/3  ea  ■ 

dXj 

l£ 

1  a0 

1 

-  eijkfjuk  and 

(Ala) 

do 

d 

dt  ~ 

--(UjO  +  U 

(Alb) 

with  the  velocity  subject  to  the  incompressibility  con¬ 
straint 
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In  (Al)  and  (A2)  spatially  filtered  variables  are  de¬ 
noted  by  (  ).  Other  variables  in  (Al)  include  the  rota¬ 
tion  vector  f  =  (0,  0,  /)  where  /  is  the  Coriolis  param¬ 
eter,  gravity  is  g,  the  reference  temperature  and  density 
are  ( 0o ,  p„),  and  the  generalized  pressure  p*  =  plpa  + 
(%)e,  where  e  =  tuI2  is  the  subgrid-scale  energy.  The 
large-scale  externally  imposed  pressure  gradients 


~7M  =  (-fv-fu-0)  (A3) 

are  prescribed  in  terms  of  the  x-y  components  of  the 
geostrophic  wind  ( Ug ,  Vg).  The  pressure  p*  at  each  time 
step  is  the  solution  of  a  Poisson  equation  (Sullivan  et  al. 
1996)  formed  from  the  discretized  continuity  Eq.  (A2). 

As  a  consequence  of  spatial  filtering,  subgrid-scale 
(SGS)  fluxes  of  momentum  t y  =  ujiij  —  TiiTij  —  (2/3)Sye 
and  its  scalar  i j/,-  =  —  0z7,  appear  in  (Al).  In  the 

present  LES  model  these  unknown  fluxes  are  modeled 
using  simple  eddy  viscosity  prescriptions 
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with  eddy  viscosity  v,  =  ckley2\  ck  is  a  modeling  constant 
and  the  length  scale  /  is  set  equal  to  the  LES  filter  A, 
except  in  regions  of  stable  stratification  where  it  is  re¬ 
duced  (Deardorff  1980).  We  note  that  SGS  modeling  is 
an  active  research  area  and  numerous  alternate  ap¬ 
proaches  to  modeling  these  SGS  variables  are  available 
(e.g.,  Meneveau  and  Katz  2000;  Geurts  2001;  Sullivan  et 
al.  2003;  Wyngaard  2004;  Sullivan  et  al.  2006a;  Hatlce 
and  Wyngaard  2007).  LES  solutions  near  flat  param¬ 
eterized  z0  boundaries  are  dependent  on  the  SGS  clo¬ 
sure  (e.g.,  Mason  and  Thomson  1992;  Sullivan  et  al. 
1994;  and  others).  This  dependence  is  expected  to  be 
weaker  in  the  present  application  as  the  surface  form 
(pressure)  stress  is  resolved.  In  other  words  there  is  less 
reliance  on  the  SGS  model  to  support  surface  fluxes  in 
the  presence  of  resolved  wavy  surfaces.  This  specula- 
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lion  is  supported  by  direct  numerical  simulations  Sulli¬ 
van  et  al.  (2000),  but  clearly  requires  further  research. 

The  transport  equation  for  subgrid-scale  energy  e  = 
t,,/2  is  (Deardorff  1980) 


de  d  _  T jj  /  dut  dUj\ 

dt  dXj^Uie^  2  \dXj  +  dxj 
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In  anticipation  of  the  numerical  solution  of  (Al)  and 
(A2)  we  introduce  contravariant  flux  velocities 
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which  point  in  directions  perpendicular  to  the  compu¬ 
tational  cell  faces  (tj-£,  £-rj),  respectively.  They 

satisfy  the  transformed  continuity  equation 


and  includes,  in  order,  time  tendency,  advection,  SGS 
production,  buoyancy,  diffusion,  and  viscous  dissipa¬ 
tion.  The  latter  is  modeled  using 

g3/2 

®  =  C‘D~\~  ■  (A6) 

The  subgrid-scale  constants  ( ck ,  c®)  (0.1,  0.93) 

(Moeng  and  Wyngaard  1988)  are  determined  by  apply¬ 
ing  a  sharp  cutoff  filter  in  the  inertial  subrange  and 
matching  with  a  Kolmogorov  spectrum. 

The  specific  choice  of  grid  transformation  from 
physical  to  computational  space  relies  on  the  problem 
definition.  We  assume  the  underlying  wavy  surface  is 
an  externally  imposed,  two-dimensional,  single  plane 
wave  of  height  h(x,  t)  =  a  cos[k(x  —  ct)]  with  amplitude 
a.  wavelength  A  (or  wavenumber  k  =  2ttIX)  propagating 
with  phase  speed  c.  The  wave  is  further  assumed  to 
obey  the  linear  dispersion  relationship  c2  =  gXI2ir.  As 
in  second-order  closure  modeling  (e.g.,  Gent  and  Tay¬ 
lor  1976)  we  introduce  a  streamwise  coordinate  x'  = 
x  —  ct  to  freeze  the  movement  of  the  surface.  The  time 
and  streamwise  advective  operators  in  (Al)  for  any 
field  g  then  transform  as 
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Allowing  the  grid  to  advect  horizontally  with  transla¬ 
tion  speed  c  is  equivalent  to  applying  a  Galilean  trans¬ 
formation  to  the  governing  equations:  the  model  equa¬ 
tions  are  not  fundamentally  changed  if  we  replace  Ti  by 
u'  =  u  —  c.  Next,  the  physical  space  coordinates  (x\  y, 
z )  are  mapped  to  computational  coordinates  £  =  = 

(£,  t),  4)  using  the  grid  transformation 

£=£(x',z),  r 7=y,  £=£(x',z),  (A8) 


with  the  Jacobian  of  the  transformation  /  =  —  £z£x. 

In  transforming  the  flow  equations,  described  below, 
we  frequently  make  use  of  the  grid  transformation 
identity  (Anderson  et  al.  1984) 
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Applying  the  grid  transformation  (A8)  to  (Al)  and 
(A5),  and  making  use  of  the  continuity  Eq.  (All)  and 
the  identity  (A9)  yields  the  transformed  set  of  LES 
equations  in  strong  conservation  form: 
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b.  Boundary  conditions 

The  flow  is  assumed  to  be  homogeneous  in  horizon¬ 
tal  planes  and  thus  explicit  boundary  conditions  only 
need  to  be  specified  at  the  water  surface  and  at  the  top 
of  the  computational  domain.  The  upper  boundary  in 
physical  space  is  located  far  above  the  wavy  surface  and 
as  a  result  the  horizontal  gridlines  are  flat  and  orthogo¬ 
nal  to  vertical  grid  lines.  At  the  top  of  the  computa¬ 
tional  domain  the  upper  boundary  conditions  are  sim¬ 
ply  set  equal  to  the  values  in  our  flat  LES  code.  We  use 
a  radiation  condition  (Klemp  and  Durran  1983)  for  ver¬ 
tical  velocity  along  with  a  specified  constant  gradient 
for  potential  temperature,  zero  vertical  gradients  for 
the  horizontal  velocities,  and  zero  SGS  turbulence 
fields  (Moeng  1984). 

The  important  changes  occur  at  the  lower  boundary 
z  =  Zb  =  h(x,  t)  where  the  Cartesian  velocity  compo¬ 
nents  are  set  equal  to  the  orbital  velocity  of  the  re¬ 
solved  wave  («,  iv)x  =  akc{cos[k(x  —  cf)],  sin[k(x  —  ct)]) 


=  0  for  i  =  1,  2,  3. 
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[e.g.,  section  3.2  of  Lighthill  (1978)].  In  the  frame  of 
reference  moving  with  speed  c  the  boundary  conditions 
on  the  Cartesian  velocity  components  are 

-c  +  u~ 

0  ,  (A13) 

w 

rvs  _l 

which  requires  the  contravariant  flux  velocities  satisfy 

~  JU~  I"  (~c  +  us)£x  +  wsCx~ 

JV  =  0  .  (A14) 

_JW  \  |_(-c  +  us)$z  +  wsCz_ 

For  a  sinusoidal  waveform  of  small  wave  slope  ak  C  1 
the  boundary  conditions  on  the  contravariant  flux  ve¬ 
locities  become  /  ([/,  V,  W)  — »  (— c,  0,  0)  (Sullivan  et  al. 
2000). 

A  key  difference  between  the  LES  and  DNS  is  the 
formulation  of  the  boundary  condition  on  surface 
fluxes.  The  LES  is  intended  to  model  a  high  Reynolds 
number  geophysical  flow  and  in  this  regime  it  is  not 
computationally  feasible  to  resolve  the  viscous  sublay¬ 
er.  Therefore  we  apply  bulk  drag  formulas  to  estimate 
subgrid-scale  surface  fluxes  of  momentum  and  scalars, 
similar  to  the  approach  in  an  LES  with  a  flat  boundary. 
The  total  drag  on  the  PBL  consists  of  resolved  form 
stress  and  an  unresolved  viscous  drag  that  rides  on  the 
wavy  surface.  In  the  LES,  a  high-Reynolds  number  sur¬ 
face  drag  law  based  on  a  “ zQ ”  boundary  condition  is 
adopted  at  the  lower  boundary,  essentially  a  law-of-the- 
wall  expression  is  applied  instantaneously  at  every  sur¬ 
face  grid  point  to  relate  the  surface  winds  and  fluxes. 
Mason  and  Callen  (1986)  and  Wyngaard  et  al.  (1998) 
discuss  the  applicability  of  this  approximation,  which  is 
experimentally  verified  in  a  flat  plate  boundary  layer 
flow  by  Nakayama  et  al.  (2004).  We  adapt  this  law-of- 
the-wall  parameterization  to  our  wavy-surface  applica¬ 
tion,  similar  to  the  methodology  used  in  second-order 
closure  modeling  (Gent  and  Taylor  1976;  Li  1995).  In  a 
neutrally  stratified  flow,  the  surface  friction  velocity  uA. 
due  to  the  unresolved  surface  waves  (or  roughness)  is 
estimated  from 

u*  I  A£/2\ 

111^/2)1  =  -ln(^— j,  (A15) 

where  k  =  0.4  is  the  von  Karman  constant,  za  is  the 
specified  surface  roughness,  and  A£  =  ^  —  h  is  the 
normal  distance  from  the  surface  wave  h  to  the  first  £ 
grid  point;  Uy  =  (u  ■  s),  where  s  and  lUyl  are  the  wind 
vector  and  wind  speed  parallel  to  the  surface,  with  s  the 
unit  vector  tangent  to  the  surface.  The  surface  momen¬ 
tum  flux  r  estimated  from  the  bulk  formula 


assumes  the  surface  stress  and  the  surface  wind  are 
parallel.  In  PBL  flows  with  surface  heating  or  cooling  a 
constant  surface  buoyancy  flux  is  specified,  and  u *  is 
modified  using  Monin-Obukhov  similarity  functions 
(Moeng  1984).  This  correction  is  small  in  the  present 
application  since  the  first  vertical  grid  level  is  quite 
close  to  the  surface,  A£,  «=  1  m.  The  wave-following 
surface  stresses  t  are  the  physical  components  of  a  sec¬ 
ond-order  tensor  and  are  converted  into  the  Cartesian 
stress  t ij,  needed  in  (A12),  using  standard  transforma¬ 
tion  rules  (see  section  13.3  of  Wylie  1966). 


c.  Numerical  method  and  grid  generation 

The  numerical  algorithm  used  to  integrate  the  LES 
model  Eqs.  (All)  and  (A12)  is  identical  to  that  in  our 
DNS  code  (Sullivan  and  McWilliams  2002;  Sullivan  et 
al.  2000).  Lor  our  mixed  finite-difference  pseudospec- 
tral  differencing  scheme  a  special  arrangement  of  vari¬ 
ables  is  employed.  The  Cartesian  velocity  and  scalar 
variables  (u,  0,  p*,  e)  are  colocated  at  cell  centers  while 
the  contravariant  flux  velocities  (U,  V)  are  located  at 
cell  centers  with  W  located  at  cell  faces.  The  positioning 
of  Ut  mimics  the  arrangement  of  variables  in  our  flat 
Cartesian  LES  code.  Advantages  of  the  colocated  grid 
structure  are  as  follows:  1)  all  advective  terms  can  be 
compactly  discretized  using  a  skew  symmetric  form, 
namely,  [d([/;w,)/3£;  +  (7;3m;/3£;]/2  for  momentum  ad- 
vection  and  [3([/7-0)/0£,  +  t/,50/3^/2  for  scalar  advec- 
tion;  and  2)  the  location  and  orientation  of  U  maintains 
tight  velocity-pressure  coupling  as  the  continuity  equa¬ 
tion  dUfld^i  is  used  to  construct  the  discrete  pressure 
Poisson  equation.  The  spatial  discretization  is  pseu- 
dospectral  along  lines  of  constant  £  or  tj  and  second- 
order  finite  difference  in  the  vertical  coordinate  C  A 
third-order  Runge-Kutta  time-stepping  scheme  operat¬ 
ing  with  a  fixed  CFL  number  is  employed  (Sullivan  et 
al.  1996).  An  important  difference  from  a  flat  Cartesian 
code  is  the  appearance  of  variable  coefficients  in  the 
pressure  Poisson  equation.  This  prevents  a  direct  solu¬ 
tion  using  Fourier  transforms  and  tridiagonal  matrix 
inversion  for  each  pair  of  horizontal  wavenumbers. 
Here  we  use  an  iterative  solution  method  for  the  pres¬ 
sure  described  in  Sullivan  et  al.  (2000).  The  entire  code 
is  parallelized  using  the  Message  Passing  Interface 
(MPI)  with  domain  decomposition  in  £.  A  custom-built 
MPI  matrix  transpose  is  used  in  the  solution  of  the 
pressure  Poisson  equation. 

The  final  element  in  our  computational  procedure  is 
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the  generation  of  an  acceptable  field  grid.  Since  the 
underlying  waveform  is  simple  and  stationary  in  com¬ 
putational  space  an  adequate  mesh  can  be  created  using 
conformal  techniques.  Given  periodicity  in  the  horizon¬ 
tal  direction,  a  flat  upper  boundary,  and  a  specified 
surface  wave,  we  then  solve  two  standard  elliptic  grid 
generation  equations  for  the  (x',  z)  coordinates 
(Thompson  et  al.  1985;  Dimitropoulos  et  al.  1998): 


dV 


dV 


d2z 


d2z 


—  +  —  =  0  and  — +  — =  0.  (A17) 


d? 


df  K 


The  vertical  boundary  conditions  at  the  surface  £  =  0 
and  at  the  top  of  the  computational  domain  £  =  zL  are 

dx'  dh  dz 

z  =  h(x'),  -^=-^  at  £=0  and 

(A18a) 

dx' 

z  =  Zl,  =  0  at  £  =  zL.  (A18b) 


The  numerical  solution  of  these  elliptic  equations  gen¬ 
erates  smoothly  varying  grids.  The  grid  metrics  and  Ja¬ 
cobian  are  constructed  numerically  from  the  one-to- 
one  mapping  between  (x',  y,  z)  and  (£,  tj,  £).  We  note 
the  use  of  a  conformal  grid  is  quite  advantageous  in 
DNS  as  it  greatly  streamlines  the  viscous  term,  but  does 
not  lead  to  the  same  simplification  in  LES  since  the 
subgrid  flux  terms  contain  spatially  varying  eddy  vis¬ 
cosity  and  diffusivity.  Last,  in  order  to  focus  the  grid 
near  the  surface  the  vertical  spacing  is  varied  using  con¬ 
stant  algebraic  stretching;  that  is,  the  ratio  of  any  two 
adjacent  vertical  cells  is  held  constant,  K  =  A£j+1/A£;. 
Stretching  factors  K  <  1.036  vary  the  grid  smoothly  but 
at  the  same  time  provide  adequate  leverage  to  span  a 
large  vertical  extent. 
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